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Abstract 



u 

Motifs are patterns of subgraphs of complex networks. We studied the impact of 
such patterns of connectivity on the level of correlated, or synchronized, spiking ac- 
tivity among pairs of cells in a recurrent network model of integrate and fire neurons. 
For a range of network architectures, we find that the pairwise correlation coefficients, 
averaged across the network, can be closely approximated using only three statistics 
J> of network connectivity. These are the overall network connection probability and the 

frequencies of two second-order motifs: diverging motifs, in which one cell provides 
input to two others, and chain motifs, in which two cells are connected via a third 
intermediary cell. Specifically, the prevalence of diverging and chain motifs tends to 
increase correlation. Our method is based on linear response theory, which enables us 
to express spiking statistics using linear algebra, and a resumming technique, which ex- 
trapolates from second order motifs to predict the overall effect of coupling on network 
• • correlation. Our motif-based results seek to isolate the effect of network architecture 

perturbatively from a known network state. 
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1 Introduction 



Neural networks are highly interconnected: a typical neuron in mammalian cortex receives 
on the order of a thousand inputs. The resulting collective spiking activity is characterized 
by correlated firing of different cells. Such correlations in spiking activity are the focus of a 
great deal of theoretical and experimental work. This interest arises because correlations can 
strongly impact the neural coding of information, by introducing redundancy among different 
(noisy) neurons, allowing noise cancellation effects, or even serving as additional "channels" 



of information 


Zohary et al. , 1994 Singer and Gray 


1995, Shadlen and Newsome 


1998b 


Panzeri et al. 


1999 


|Abbott and Dayanl |1999l |Sompoli 


nsky et al.[ |2001[ |Panzeri et al 


|2001 


Schneidman et al. 


2003, Latham and Nirenberg, 2005 


Josic et al. , 2009 Beck et al. 


2011 . 


Correlations can also help gate the transfer of signals from one brain area to another 


Salinas 



, 2000, 


Fries 


2005 


Bruno 



1], and determine the statistical vocabulary of 
a network in terms of the likelihood that it will produce a particular spike patterns in its 
repertoire 



Schneidman et al. , 2006 , Shlens et al. , 2006 



Thus, the possible effects of correlations on the encoding and transmission of signals are 
diverse. Making the situation richer still, correlations can depend on characteristics of the 
"input" signals (or stimuli) themselves. The result is an intriguing situation in which there 
is no simple set of rules that determines the role of correlations - beneficial, detrimental, or 
neutral - in neural computation. Moreover, despite major progress, the pattern and strength 



of correlations in neuronal networks in vivo remain under debate Ecker et al. 2010, Cohen 



and Kohn, 2011 . This demands tools that will let us predict and understand correlations 
and their impact in specific types of neural circuits. 

In this paper, we take one step toward this goal: We determine how simple statistics of 
network connectivity contribute to the average level of correlations in a neural population. 
While other factors such as patterns of correlations in upstream neural populations and the 
dynamical properties of single cells contribute in important ways, we seek to isolate the role 
of connection statistics in the most direct way possible. We return to the question of how 
our results on connectivity might be combined with these other factors at several places in 
the text, and in the Discussion. 

We define the statistics of network connectivity via motifs, or subgraphs that are the 
building blocks of complex networks. While there are many ways to characterize connec- 
tivity, we choose the motif-based approach for two main reasons. First, we wish to follow 
earlier work in theoretical neuroscience in which the frequency of small motifs is used to 
define low-order statistics of network adjacency matrices Zhao et al. , 2011 . Second, recent 
experiments have characterized the frequencies with which different motifs occur in biologi- 
cal neural networks, and found intriguing deviations from what we would expect in the case 



of unstructured, random connectivity Song et al. , 2005 , Haeusler et al. , 2009 



Fig. [T] depicts the motifs we will consider: a single cell projecting to two downstream cells 
(diverging motif), a single cell receiving inputs from two upstream cells (converging motif), 
and a series of two connections in and out of a central cell (chain motif). To assess how 
prevalent these motifs are in a given network, we count their occurrences, and compare the 
observed motif counts with those expected in a random graph Song et al. , 2005 . This is 
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a regular network which has the same total number of connections, and in which each cell 
has the same, evenly divided number of incoming and outgoing connections - i.e, the same 
in and out degree. (Importantly, the (relative) motif counts for such regular graphs agree 
with those in the classical model of random graph, the Erdos-Renyi model in the limit of 
large network size; thus, when we refer to the prevalence of network motifs, this means in 
comparison to either a regular or Erdos-Renyi graph.) 




Figure 1: (Left) Counting motifs in a network. Example of a diverging motif in the network 
is shown in bold solid line. Similarly a converging and a chain motif are shown in dashed 
line and dash-dotted lines, respectively. In total, the network has 8 connections, 6 diverging, 
7 converging and 5 chain motifs. (Right) The different types of second order motifs. 



Fig. [2] illustrates the importance of network motifs in determining the average correlation 
across the network. Here, we simulate 265 different networks of excitatory and inhibitory, 



exponential integrate and fire cells Dayan and Abbot, 2001 . Importantly, we bias each 



neuron so that it fires with the same rate, regardless of its connectivity. We explain in more 
detail below that this is important to isolate the effects of network connectivity alone. The 
black dots show the average correlation for Erdos-Renyi networks that have different connec- 
tion probabilities p. As expected, correlations increase with connection probability. Next, 
the grey dots show the average correlation in networks that all have the same connection 
probability (p = 0.2), but with a different prevalence of motifs compared to the correspond- 
ing Erdos-Renyi model. Interestingly, the range of correlation values obtained at this fixed 
connection probability p is as large as that obtained in the Erdos-Renyi networks over a 
range of p values that reaches to roughly twice the connectivity. Thus, motifs - over and 
above net connectivity - play a strong role in determining the strength of correlations across 
a network. 
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Mean Correlation 
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Figure 2: Impact of changing motif frequencies on mean correlation in non-Erdos-Renyi 
networks (gray dots) compared with the effect of changing connection probability in Erdos- 
Renyi networks (black dots). Mean correlation coefficients (averaged over all cell pairs in the 
network) are plotted against connection probability p stat in the Erdos-Renyi model. Keeping 
Pstat = 0.2, and varying second order motif frequencies strongly affects average correlations. 
The curve presents the theoretical predictions from Eq. (45). For network parameters see 
Fig.|3} 



But which motifs contribute? And are the second-order motifs of Fig. [T] sufficient, or 
must higher-order motifs (involving four or more cells) be included to understand average 
correlations? Fig. [3] suggests part of the answer. Here we used the same 265 network samples 
represented by the gray dots in Fig. |2j Connection probability was kept at p ~ 0.2, but 
frequencies of second order motifs differed. First, panel A of Fig. [3] shows the dependence 
of average network correlations on the total motif counts in the network. Disappointingly, 
no clear relationship is evident. 

A more careful approach is presented in panel B of Fig. [3j First, we plot only the mean 
correlation among cell pairs of a given type - here, excitatory cell pairs. Second, we do not 
lump together all motifs shown in Fig. [1} Instead, we calculate weighted motif counts that 
are classified according to the types of constituent neurons. The three panels now exhibit 
a much clearer trend: Correlation levels vary systematically with weighted motif counts for 
the diverging and chain motifs, while there is no clear dependence on the converging motif. 
This improved linear fit is quantified by the coefficient of determination R 2 (as usual, the 
squared correlation coefficient), for the mean correlation among excitatory cells regressed 
against the motif counts. 

Our goal in the balance of this paper is to explain two aspects of the relationship suggested 
in panel B of Fig. |3j First, we show how to derive analytically a relationship between motif 
counts and network-wide correlation that successfully identifies these trends. We previously 

to derive an explicit expression for the pairwise 



built on the work of Lindner et al. 2005 



correlation between cells in a network Trousdale et al., 2012 ; close analogs of this expression 



have been derived via related approaches for Hawkes processes Hawkes, 1971b , Pernice et al. 
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2011, 2012 . This expression shows how patterns of network connectivity, together with the 
dynamical properties of single cells, shape correlations. Here we extend our previous work 
to show how the same expression can reveal the impact of motifs on correlations. 

Secondly, we use this analytical approach to obtain a result that is not readily apparent 
from Fig. |3} We show that - for a range of different network models (see Sec. [3] ) - the 
average network correlation can be closely approximated using the connection probability 
and frequencies of motifs that involve only two connections between cells (See Fig. [TJ. Specif- 
ically, we show that the prevalence of diverging and chain motifs tends to increase average 
correlation (in excitatory only networks), while converging motifs have no effect (in either 
single population or multi-population networks) . 

The paper is structured as follows. We explain in Sec. [2] our setup and introduce the 
linear response theory which is the basis of our analysis. Sec. [3] defines the three types of 
motifs and their frequencies, and discusses the classes of networks to which we test and apply 
our theory. For simplicity, we first demonstrate our analysis in a case of a population of a 
single population of excitatory cells in Sec.|4j and then generalized to interacting populations 
of excitatory and inhibitory cells in Sec. [5] We conclude by discussing a crucial assumption 
of our network model (Sec. [6]), and compare our theory with simulations of integrate and 
fire neuron networks (see Sec. [7]). Biological interpretations, limitations, and extensions are 
covered in our Discussion, and an appendix and table of notation follow. 
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Mean Regression 

Chain (x10 3 ) Chain (x 10 3 ) Converging (x 10") Correlation Coefficients 




Figure 3: Mean correlation as a function of motif frequency in networks of excitatory and 
inhibitory neurons. A Correlation coefficient averaged over all cell pairs as a function of 
total motif counts. B Correlation coefficients averaged over pairs of two excitatory cells, 
as a function of weighted sums of subgroup motif counts (see Fig. [i] and Eq. p4|) . The 
linear combination coefficients are determined by the resumming theory (see Sec. |5.3[ ). The 
network consists of 51 excitatory neurons and 49 inhibitory neurons with excitatory and 
inhibitory coupling strengths set at 22.82 mV-ms or -22.82 mV-ms. All neurons have the same 
uncoupled cellular dynamic parameters (see Eq. Q): Tj = 20,v t h = 20, v r = — 60,r re j = 
2,E Lti = -60, Ei = 8, of = 12, = -53, A T = 1.4, = 5,r Ai = 10. Spike count 
correlation coefficients were calculated using a 500ms window size. These graphs are re- 
sampled from 4096 degree distribution generated graphs using Latin hypercube sampling 
based on parameter space defined by axis in panel B over a reliably sampled region. In 
panel A, some of the motifs are given a "positive" count - i.e., if they involve two excitatory 
or two inhibitory connections in a chain - and others are given a negative count. For example 
a converging motif with one excitatory presynaptic cell and one inhibitory presynaptic cell 
will contribute negatively. 

2 Neuron models, cross-spectra, and measures of col- 
lective network activity 

In this section we describe a model spiking network composed of integrate and fire (IF) 
neurons. We introduce the measures that we will use to quantify the level of correlations 
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between cell pairs in the network. The analytical approximations of these correlations given 
in Eq. ^ will be the basis for the subsequent analysis. While we develop the theory for 
networks of IF neurons, the main ideas are applicable to many other models (e.g. Hawkes 
model). We return to this point in the Discussion. 

2.1 Networks of integrate-and-fire neurons 

In a network of integrate-and-fire (IF) units, the dynamics of each cell are described by a 
single membrane voltage variable Vi, which satisfies 



TiVi 



-(vi - E Lji ) + + Ei + Jafr&tb) + hit) . 



Here, Elj is the leak reversal potential, ip{vi) is a spike generating current, and Ei is the 
mean external synaptic input from sources not modeled. In numerical simulations we use the 
exponential integrate-and-fire (EIF) model Fourcaud-Trocme et al. 2003 , so that ip(v) = 
A T exp[(t> — Vt)/A t ]. Each cell has independent fluctuations due to internal noise and 
external inputs (e.g., from a surrounding network that is not explicitly modeled). We d escribe 
such effects by additive terms, y/ofTi£i(t), which are Gaussian white noise processes, 



White 



et al. , 2000 , Renart et al. , 2004 . Synaptic input to cell i from other cells in the network is 
denoted by fi(t) (see below). 

When the membrane potential reaches a threshold, Vth, an action potential, or spike, is 
generated; the membrane potential is then reset to a lower voltage v r and held there for 
an absolute refractory period r r . We denote the time at which the j th neuron fires its k th 
spike as t^; taken together, these define this neuron's spike train yj(t) = ^2 k S(t ~ tj,k) ■ 
Importantly, synaptic interactions among neurons are initiated at these spike times, so that 
the total synaptic input to the i th cell is 



/*(*) = 5^( J w * where J ii(*) 







exp 



t—TD, 



t > T DJ 
t < T D ,j 



In the absence of a synaptic connection from cell j to cell i, we set W,y = 0. Hence, the 
N x N matrix of synaptic weights, W, defines a directed, weighted network. 

In simulations we used the parameters given in the caption of Fig. [3j unless stated 
otherwise. 



2.2 Measure of network correlation 

Dependencies between the responses of cells in the network may be quantified using the spike 



train auto- and cross-correlation functions Gabbiani and Cox, 2010 . For a pair of spike 



trains yi(t),yj(t), the cross-covariance function is given by 

C ij (r) = E[(y l (t)-r i )(y j (t + r) 



rj)) 
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where and Tj are the firing rates of the two cells. Here and throughout the paper we 
assume that the spike trains form a (multivariate) stationary process. The auto-correlation 
function is the cross-correlation of the output of the cell and itself, and C(t) is the matrix 
of cross-correlation functions. The matrix C(u) with entries defined by 



is the matrix of cross-spectra (see next section). The cross-spectrum of a pair of cells is 



equivalent to the Fourier transform of their cross-correlation function. Stratonovich , 1967 
Denote by N yi (ti,t 2 ) = f* 2 yi(s)ds the spike count of cell % over a time window [ti,t 2 ] 
The spike count correlation p^iT) over windows of length T is defined as 



cov (N yi (t,t + T),N yj (t,t + T)) 
var (N m (t, t + T)) var (N y . (t, t + T)) 



We will make use of the total correlation coefficient py(oo) = limj-^oo Pij{T) which captures 
dependencies between the processes y^Uj over arbitrarily long timescales, but may also 



describe well the nature of dependencies over reasonably short timescales Bair et al. , 2001 



de la Rocha et al. 2007, Shea-Brown et al. 2008, Rosenbaum and Josic 2011 . The spike 



count covariance is related directly to the cross-correlation function by Tetzlaff et al., 2008 



cov (N yi {t,t + T),N y .{t,t + T)) = [ Cij(r)(T — \r\)dr. 

J-T 

Thus, total correlation may be defined alternatively in terms of the integrated cross-correlations 
(equivalently, the cross-spectra evaluated at co = 0): 



C<,(0) 



C«(0)Ctf(0) 



(2) 



Throughout this paper, we will use p i ,(oo) as the measure of correlation and use the above 



equation to calculate it from the cross-spectra matrix Pernice et al. 2011 . However, our 



analysis can be similarly applied to study Cy(u;), and hence the entire correlation function 



in time Trousdale et al., 2012 . 



2.3 Linear response approximation of cell response covariance 

Linear response theory 



Gabbiani and Cox, 2010, Risken, 1996 can be used to approxi 



mate the response of sin 


gle cells, and the joint response of cells in a network |Lindner and 


Schimansky- Geier 


2001 


Brunei et al. , 


2001 


Lindner et al. , 


2005 


Trousdale et al. 


2012 . 



Consider an IF neuron obeying Eq. (|lj), but with the mean of the inputs f(t) absorbed into 
the constant E. We denote the remaining, zero-mean input by eX(t), so that 

rv = -(v - E L ) + ij(v) + E + v^Ve(^) + eX(t) . (3) 
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For fixed input fluctuations eX(t), the output spike train will be different for each realization 
of the noise and each initial condition v (0). The time-dependent firing rate is obtained 
by averaging the resulting spike train over noise realizations and a stationary distribution 
of initial conditions. For all values of e, this stationary distribution is taken to be the one 
obtained when e = 0. We denote the resulting averaged firing rate as r(t) = (y(t)). Linear 
response theory approximates this firing rate as 

r(t)=r + (A*eX)(t), 

where r is the firing rate in the absence of input (e = 0), and the linear response kernel, A(t), 
characterizes the response to first order in e. This approximation is remarkably accurate over 



a wide range of parameters; for example, see Ostojic et al. , 2009, Richardson, 2009 



Next, we turn to the problem of approximating the output of a cell on a single trial, 
rather than the average across trials. We denote the Fourier transform of a function / 
by / = J-(f). However, for spike trains we adopt a convention that tjiiu) is the Fourier 



transform of the mean subtracted spike train yi(t) — jy Following Lindner and Schimansky 



Geier, 2001, Lindner et al. , 2005 Trousdale et al. 2012 we approximate the spiking output 



of a cell self-consistently by 



yi\w) ~ y 




(4) 



where y^{oS) is a realization of the output of cell i in the absence of input. Defining the 
interaction matrix K with entries K^(t) = {Ai * Jy)(t), we can use Eq. Q to solve for the 
vector of Fourier transformed spike train approximations 



y(w) = (I - K(w))"Y(w) 



and matrix of cross-spectra 



C(w) « C°» = (I - K(w))- 1 (»°(a;)y *(a;)>(I - K*(w)) _1 
= (I - AWF) _1 C°(I - F*W T A*)- 1 • 



(5) 



(6) 



Here A, C°, F are diagonal matrices: Aa(u) = Ai(u) is the linear response of cell i, C° (u) = 
Cf(u) = (y° ((jj)Vi* (co)) is its "unperturbed" (i.e., without coupling) power spectrum, and 
Fjj(o;) = Fi(ui) is the Fourier transform of the synaptic coupling kernel from cell i. As noted 
later, our results and analysis will hold at all frequencies and thus can be used to study 
correlations at all timescales. The weighted connectivity matrix W, defines the structure of 
the network. 

To simplify the exposition, we initially assume certain symmetries in the network. For 
instance, we consider homogeneous networks in which cells have identical (unperturbed) 
power spectra, linear response functions, and synaptic kernels. In this case the diagonal 
matrices in Eq. ^ act like scalars. We slightly abuse notation in this case, and replace 
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Ai(u) by A(cj), and C^ipj) by C°(u). This allows us to disentangle the effects of network 
structure from the effects of neuronal responses on network activity. The resulting cross- 
spectrum matrix (evaluated at u = 0) is 

C 0O = C°(l-iw)~ 1 (I-iW T )" 1 (7) 

(note that i^(0) = 1 by definition). We use this simpler expression in what follows, and 
return to heterogeneous networks in Sec [6] Since we consider only total correlation, we omit 
the dependence of the spike count correlation coefficient on window size T. Additionally, all 
spectral quantities are evaluated at w = 0, so we also suppress the dependence on u>. Finally, 
we define average network correlation by 

i N 

In subsequent sections, we will examine the average covariance across the network 

1 - 

This average cannot be directly related to that in Eq. (pi), where individual summands are 
normalized, and diagonal terms are excluded. The motif-based theory we develop predicts 
(C°°), and gives no information about the specific entries C°° . However, p avg can be deter- 
mined approximately from (C°°) alone. We describe these approximations in Appendix |A| 



2.4 Applicability of linear response theory 

Our methods depend on two, related sets of conditions for their validity. First, we take our 
cells to be driven by a white noise background. This background linearizes the response of 
the cells to sufficiently weak perturbations, improving the accuracy of the approximation 
Q; its presence is our first condition. 

Second, turning to network effects, we assume that the spectral radius ^(K) < 1, which 
gives non-singularity of the approximating processes in Eq. Q and allows us to make the 



series expansion we describe in Sec 4.2 In practice, we have found that the linear approx- 
imation to correlations will cease to provide an accurate approximation before this occurs, 
likely owing in part to a failure of the perturbative approximation. Furthermore, the ap- 
proximation seems to be most accurate at weak interaction strengths as characterized by a 
small radius of the bulk spectrum of K. 



For Erdos-Renyi networks, Rajan and Abbott 2006 derived an asymptotic (large N) 



characterization of the spectral radius of the synaptic weight matrix. In particular, for an 
Erdos-Renyi network consisting of only excitatory cells with synaptic weight w, there will be 
a single eigenvalue at pNw with the remaining eigenvalues distributed uniformly in a circle 
around the origin of radius w^/p(l — p)N. In networks with both excitatory and inhibitory 
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populations, there is a single outlier at pN E WE + pNjWj and all other eigenvalues will be 
distributed (non-uniformly) within the circle of radius \/p(l — p)(New e + Njwj). We will 
use these expressions for the bulk spectrum to quantify the strength of interactions given by 
the asymptotic spectral radius of K = AW (referred to as the Erdos-Renyi spectral radius). 
We note that we used IF simulations to directly confirm the accuracy of the linear re- 



sponse approximation for excitatory-inhibitory networks in Fig. 14 For a complete discussion 



of the performance of the linear response theory, see Trousdale et al., 2012 



3 Graphical structure of neuronal networks 

Our main goal is to determine how the small-scale statistical structure of directed networks 
influences the collective dynamics they produce - namely, the strength of spike correlations 
in networks of model neurons. We will quantify network structure using the probability of 



connections between pairs and among triplets of cells, organized into network motifs Song 



et al., 2005, Zhao et al., 2011 . 



3.1 Network motifs and their frequencies 

A motif is a subgraph composed of a small number of cells. We classify motifs according to 
the number of edges they contain. We begin by considering directed networks composed of 
identical cells. First order motifs contain one connection and hence come in only one type 
- two cells with a one-way connection. Second order motifs contain two connections, and 
therefore involve at most three interacting cells. These motifs come in three types: diverging, 
converging and chain motifs (See Fig. [TJ |Song et ah 2005, Zhao et al. , 2011] . (Note that in 
our definition, a cell can appear twice in the triplet of cells that define a second order motif. 
For example, the chain motif in Fig. [T] is equivalent to a bidirectionally coupled pair of cells 
when i = j.) 

We will consider mainly the impact of second order motifs, over and above first order 
effects. The three motifs shown in Fig. [T] arise naturally in our analysis of correlated spiking 
activity. In particular, we will show that the frequency at which each motif occurs in the 
network can accurately predict levels of correlation across the network. 

We next introduce notation that will allow us to make these ideas precise. Let W° be 
the adjacency matrix, so that - = 1 indicates the presence of a directed connection from 
cell j to cell i, and W°j = indicates its absence. To quantify the frequency of a motif 
in a given graph, we first count the total number of times the motif occurs, and divide by 
the total number of possible occurrences in a graph of that size. For first order motifs this 
definition gives the empirical connection probability, 



P 



1,3 



The preponderance of second order motifs is measured in two stages. First, we similarly 
normalize the motif count. Second, we subtract the value expected in a reference graph. 
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The resulting expressions are: 



<?di 



E«)/iV 3 -p 2 = fe( W ° W0T )y) /N 3 -p* 

i,j,k \ i,j / 

(D w ° w %) /n 3 -p 2 , 



(9) 
(10) 

fin 



where W 0T denotes the transpose of W°. Consider the expression defining q^ v : the sum 
in the first equality simply counts the total number of connections from one cell (k) to two 
others (i and j), and divides by the total number of possible connections of this type {N 3 ). 
This can be written as matrix multiplication followed by a sum over all entries i, j, as shown. 
In each case we subtract the value p 2 , which corresponds to the frequency of the motif in a 
regular graph, as well as the asymptotic frequency in an Erdos-Renyi graph as the number of 
cells, N, diverges to infinity Indeed, for Erdos-Renyi graphs any edge is present with chance 
p, and any second-order motif requires the presence of two edges. Thus, gdiv corresponds 
to the propensity for a network to display diverging motifs, over and above expectations 
from an Erdos-Renyi or regular network. The other measures in Eqns. (|9 11) have similar 
interpretations. 

The quantities in Eqns. ([9 11) can also be interpreted as empirical measures of covariance 



Roxin, 2011, Zhao et al. 2011 . If we denote by E e [-] the empirical average over all entries 



in a given network, then we may write 

g div = E e [W» fc wy - E e [WjJ E e [W° fc ] • 

Equality is attainable for g div and q con (but not simultaneously). We also note that the 
quantities defined in Eqns. (|9 11 ) are not independent (but do have three degrees of freedom). 
If we first sum over indices i, j in Eq. (JoJ) , we can rewrite the expression in terms of in and 
out degrees, {d 1 ^}, {c^ ut }. For example, 



g div = $^K ut ) 7^ 3 -P 2 = ™r(d out )/N 2 



is the scaled sample variance of the out degree across the network. Similarly, 

q COD = var(cT)/iV 2 , q ch = cov(d out , cT)/iV 2 , 



where cov(-, •) denotes sample covariance |Zhao et~aL 2011 . Therefore, 

<?div > 0, q con > 0, |g c h| < y/ <?div<?con • 

We further show in Appendix [B] that 

|?div|, |<?con|, |<?ch| < P(l -P)- 



(12) 

(13) 
(14) 
(15) 
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Eqns. (14 15) identify the attainable ranges of motif frequencies. In generating our networks, 



we compare the extent of motif frequencies that we produce using a particular scheme against 



this maximum possible range (see also Sec. 3.2) 



In networks composed of excitatory and inhibitory cells, we can represent interactions 
between cells using a signed connectivity matrix. Edges emanating from inhibitory neurons 
are represented by negative entries and those from excitatory neurons by positive entries. 
In this case, motifs are further subdivided according to their constituent cells. For instance, 
there are 6 distinct diverging motifs, since if we list all 2 3 group types of 3 cells, we see 
that the motif E <— I — > I is the same as /•<—/—>■ E, and E E — > I is the same as 
/ <r- E — >■ E. Similarly, there are 6 distinct converging motifs, and 8 distinct chain motifs, 
for a total of 20 distinct second order motifs (See Fig. El). 
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Figure 4: Second order motifs in populations of excitatory and inhibitory cells: There are 
20 subtypes of the 3 main motif types. Triangles represent excitatory neurons and circles 
for inhibitory neurons. 

This will clearly lead to some cumbersome notation! Therefore, while populations of 
interacting E and I cells are our ultimate goal, we first describe our ideas in a population of 
cells of a single type. 



3.2 Generating graphs with given motif frequency 

To numerically examine the impact of motif frequency on dynamics, we need to generate 
graphs that are equal in connection probability, but differ in the preponderance of second 
order motifs. The empirical connection probability in network samples will have small fluc- 
tuations around the statistical (i.e. expected) value we fixed. We use two ways of generating 
such graphs. The first is the degree distribution method Chung and Lu, 2002] (related to 
configuration model 



Roxin, 2011, Newman, 2003, Newman and Watts 



2001 



Here, fol- 
lowing Zhao et al.l 2011 we use a two-sided power law, with various rising and decreasing 
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exponents, peak locations and truncation limits, as expected in- and out-degree distribu- 



tions. The other is the second order network (SONET) method (for details see Zhao et al. 



2011] ) . Network samples generated using both methods cover the range of motif frequencies 
observed experimentally in cortical circuits |Song et aL 2005, Zhao et al. , 2011 . Naturally, 



this experimentally observed range is smaller than the full extent of possibly attainable fre- 
quencies (see Eqns. (14 15); however, the SONET method covers this full range as well. 



Details are given in Appendix |CJ) . 

We use both methods to generate network samples in the excitatory only case and found 
similar results; here, we only show data generated using the SONET method as it covers 
a larger range of motif frequencies. In excitatory-inhibitory networks, we use the degree 
distribution method. 

We emphasize that our approaches below do not a priori specify any particular way 
of generating network samples. However, their accuracy will depend on how this is done. 
In particular, while we find that our approach is quite accurate for the family of networks 
that we described above, they can break down in networks that have significant additional 
structure, a point to which we will return. 



4 Impact of second order motifs in networks of excita- 
tory cells 

As we have shown in Section [2j linear response theory can be used to approximate cross- 
correlations between cells in a neuronal network. We next explore how the key expression, 
given in Eq. Q, can be applied to relate the frequency of second order motifs to the average 
correlation across pairs of cells. For simplicity, we first consider networks consisting only of 
a single class of excitatory cells. The results extend naturally to the case of two interacting 
populations of excitatory and inhibitory cells, as we show in Sec. [5} 



4.1 Results: linear dependence between mean correlation and mo- 
tif frequencies 

Fig. [5] illustrates the relationship between second order motif frequency and average corre- 
lation in networks of excitatory cells. In all examples in this section, we use networks of 
N = 100 cells, with parameters as given in the caption of Fig. [3j and with graph structures 



generated by the SONET algorithm (See Sec. 3.2). Correlation coefficients between cells are 



computed using the linear response approximation given in Eq. (|6j). We find that average 
correlations depend strongly on the frequency of diverging and especially chain motifs, but 
only weakly on the frequency of converging motifs. To quantify the linear fit, we use the coef- 
ficient of determination between the linear fit and the result of Eq. ([6]), obtaining R 2 = 0.80. 
This high value suggests that the second order motifs are highly predictive of network cor- 



relation, in the simplest possible (linear) way (as explained in Appendix D, this prediction 



can be further improved when we compensate for the fluctuations in empirical connection 
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probabilities due to the finite size of network samples). In the balance of this section, we 
explain why this is the case, derive via a resumming theory a nonlinear predictor of mean 
correlation in terms of motif frequencies, and extract from this an explicit linear relationship 
between the probability of observing second order motifs and the mean correlation in the 
network. 



Mean Regression 

Chain (x 10" 2 ) Chain (x10 2 ) Converging (x 10" 2 ) Correlation Coefficients 




5 10 Converging 5 10 Diverging 5 10 Diverging a. o o 

(x10' 2 ) (x10" 2 ) (x10" 2 ) < 1 * 



Figure 5: (Left) The relationship between second order motif frequencies and average 
correlation in purely excitatory networks. Each dot corresponds to a network sample and its 
shading represents the corresponding average correlation coefficient computed using Eq. (J6|. 
Each axis represents one of the three second order motif types as defined in Eq. ([9]- 10). The 
effective coupling strength, characterized by the Erdos-Renyi spectral radius, was 0.2 for all 



networks (See Sec. 2.4). (Right) Bars show the linear regression coefficients calculated from 



the resumming theory (See Eq. (31)) and numerically (points) between motif frequencies 
and average network correlations. The error bars around each point denote 95% confidence 
intervals for the regression coefficients. The coefficient of determination R 2 is 0.80. The 
data was obtained from 512 network samples generated using the SONET algorithm (See 
Sec. 



3.2) 



4.2 First theory: second order truncation for network correlation 

If the spectral radius $(iW) < 1, then the matrix inverses in Eq. ([7]) can be expanded in 
a power series in AW as 

A i+j W\W T y (16) 



c° 



i,j=0 



Horn and Johnson, 1990 



Terms in this expansion correspond naturally to paths through the 

2009a] |b 



graph defining the network jTrousdale et al. 2012 Pernice et al. , 2011, Rangan 



For example, the terms (A 2 W 2 )ij = A 2 ^2 k WikWkj give the contributions of length two 
chains between pairs of cells Entries in the matrices A 3 W 2 W T , A 3 W(W T ) 2 give the 
contributions of diverging motifs of third order consisting of a projection from a single cell 
to a cell pair. In this case, one branch of the motif is of length two, while the other is a 
direct connection. 
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Let I7V1 7V2 denote the N% x N 2 matrix of ones, and define the N— vector L = ^1tv,i- We 
define the orthogonal projection matrices H, © which will play a crucial role in the following 
analysis: 

H = iVLL T , © = I — H. (17) 

Note that if X is an N x N matrix, then L T XL =: (X) is the empirical average of all 
entries in X. We first observe that the empirical network connection probability can be 
obtained from the adjacency matrix, W°, as 

p = L T W°L. 

We can also express second order motif frequencies in terms of intra-network averages. For 
instance, 

q div = ±L T W°W 0T L - p 2 

= Il/W° (H + @)W 0T L-p 2 

N 1 (18) 

= (L T W°L) (L T W 0T L) + ^L T W°eW 0T L - p 2 

= lL T w°ew 0T L. 

N 

Similarly, g con , g c h may be expressed as 

q con = ^L T W 0T W°L - p 2 

= I L T w or ew°L, 



(19) 



and 



g ch = ^L T W°W°L-p 2 



= 1l t w 0T w 0T l- p 2 (20) 

N 

= — L T w°ew°L. 

N 

To relate second-order motif frequencies to mean correlations between pairs of cells, we 
can truncate Eq. (16) at second order in (AW), giving 

C°° ~ ~ /~\2 / ~ \ 2 o/~\2 



i + AwW° + AwW 0T + (Aw} w°w 0T + (Aw} (w°) 2 +(Aw} (w 0T ) 2 . (21) 

To obtain the empirical average of pairwise covariances in the network, (C°°), we multiply 
both sides of Eq. (21) on the left and right by L T and L, respectively. Making use of 
Eqns. (18 20), we obtain 

(C°°) 1 ~ /~\ 2 /~\ 2 /~\ 2 

1 2Awp + 3N ( Aw) p 2 + N(Aw) q div + 2N(Aw) q ch . (22) 



C° N 
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How well does this second-order truncation predict levels of correlation across different 
neural networks? Fig. [6] shows that the truncation correctly captures general trends in levels 
of correlation from network to network, but makes a substantial systematic error. Here, 



we plot correlations predicted with the truncated Eq. (|22|) as an approximation of the full 

Indeed, the truncated 
We conclude 



expression (i.e., to all orders) for average correlations given by Eq. (J6J) 
expression gives consistent predictions only at very small coupling strengths, 
that the terms which were discarded (all terms of order three and higher in AW) can have 
an appreciable impact on average network correlation, and will next develop methods to 
capture this impact. 
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Figure 6: Scatter plot comparing the prediction of average network correlation obtained 
using Eq. ^ (horizontal axes) to the second order truncation in Eq. (22) (vertical axes). 
The diagonal line y = x is plotted for reference. Each panel corresponds to a different 
coupling strength in the same set of 512 adjacency matrices. The effective coupling strength 
is characterized by the Erdos-Renyi spectral radius shown at the top of each panel (See 



Sec. 2.4) 



4.3 Improved theory: resumming to approximate higher-order 
contributions to network correlation 

A much better approximation of the average network correlation can be obtained by con- 
sidering the impact of second order motifs on higher order terms in the expansion given by 
Eq. ( |l6| ). Note that in an Erdos-Renyi network, every motif of order m occurs with proba- 
bility p m (with the exception of motifs which involve the same connection multiple times), 
so that on average gdiv = <?con = <?ch = 0. For non-Erdos-Renyi networks, the expected values 
of gdiv, <?con and g c h are typically not zero. As we show next, the introduction of additional 
second order structure in the network also affects the frequency of motifs of higher order. 
Consider the average covariance as determined from the full order linear response ap- 



proximation of Eq. (16): 



(C°°) 
C° 



i,j=0 



Aw 



i+j 



L T (W°) l (W 0T y'L. 



(23) 
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We will first express every term in the sum given by Eq. (|23j) approximately in terms of 
first and second order motif frequencies (p, g<jiv, <?con and g c h). We illustrate this approxima- 
tion in two examples, before proceeding to the general calculation. 

Consider the term L T (W ) 3 L corresponding to the average number of length three chains 
connecting a pair of cells. Using I = H + 0, we can proceed as in the computation leading 



to Eq. (18) 



L T (W°) 3 L = L T W°(H + 0)W°(H + 0)W°L 

= L T W°HW°HW°L + L T W°eW°HW°L (24) 
+L T W°HW°eW°L + L T W°0W°0W°L. 

We next replace H by iVLL T in the last expression to obtain 

L T (W°) 3 L = iV 2 (L T W L)(L T W L)(L T W L) + iV(L T W eW L)(L T W L) 

+A^(L T W°L)(L T W°0W°L) + L T W o 0W°0W o L. (25) 

Here, the first three terms are composed of factors that correspond to the connection prob- 
ability (L T W°L) and second order chain motif frequency (L T W°0W°L). These terms 
provide an estimate of the frequency of a length three chain in a graph, in terms of the 
frequency of smaller motifs that form the chain. The last term, L T W°0W°0W o L, gives 
the frequency of occurrence of the length three chain in addition to that obtained by chance 
from second order motifs. Such higher order terms will be gathered separately and denoted 
by h.o.t. in the approximation. We therefore obtain, 

L T ( W°) 3 L = N 2 (p 3 + 2pq ch ) + h.o.t. (26) 

The exact form of this expression can be understood by referring to Fig. [7j The leading 
N 2 denotes the number of possible length three chains between a pair of cells (such a chain 
can pass through N 2 different intermediate pairs of cells) and p 3 represents the probability 
that one of these length three chains is "present" in an Erdos-Renyi graph. Recall that 
g c h measures the probability that a length two chain is present, above that expected in an 
Erdos-Renyi network. Therefore, pq c ^ represents a second order estimate of the probability 
above (Erdos-Renyi) chance that the first two connections in the chain (q c h) and the last (p) 
are present simultaneously. The prefactor 2 appears because this may also occur if the first 
and final two connections are present simultaneously. 
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Figure 7: Estimating the number of occurrences of a third order chain motif. The first 3 terms 



in Eq. (25) correspond to different ways that a chain of length three can be composed from 
smaller motifs. As explained in the text, each of the three terms represents the estimated 
probability that smaller motifs occur in one of the arrangements on the right. 

As a second example, consider the term L T (W ) 2 (W 0T ) 2 L, corresponding to indirect 
diverging motifs with two connections on each branch. A computation similar to that used 



to obtain Eq. (26) now gives 



L T (W°) 2 {W 0T ) 2 L = N 3 [p 4 + p 2 q div + 2p\ h + q 2 ch ] + h.o.t. (27) 

The four terms in this sum can be understood with the help of Fig. [8j p 4 represents the 
chance of observing the motif in an Erdos-Renyi network. The product p 2 q<nv is the estimated 
probability above (Erdos-Renyi) chance that the motif is formed by two connections ema- 
nating from the source cell, present simultaneously and independently with two connections 
emanating from the tips of the branches. The term p 2 q c h gives the estimated probability, 
above (Erdos-Renyi) chance, that one branch is present, simultaneously and independently, 
from two connections which form the other branch (the prefactor 2 concerns the probabil- 
ity of this occurring in each of the two branches). The last term, q 2 h , gives the estimated 
probability, above the Erdos-Renyi chance level, that two length two chains simultaneously 
emanate from the root cell. 



In Eq. (27), h.o.t. denotes the three distinct terms 



L T w°ew°ew 0T ew 0T L, n(l t w o l)(l t w o 0w ot 0w ot l), 
iv(L r w ew ew 0T L)(L T w 0T L). 

These terms contain two or more occurrences of in one factor, and hence correspond to 
motifs of higher than second order. In general, a factor which contains m occurrences of © 
will depend on the frequency of a motif of order m + 1, beyond that which is imposed by 
the frequency of motifs of order m. 
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Figure 8: Estimating the number of occurrence of a fourth order motif. Eq. (27) can be 



understood by decomposing this motif into the constituent first and second order motifs. 



The idea behind these two examples extends to all terms in the series in Eq. (23), as- 



suming absolute convergence. Each term in the resulting series contains a factor of the form 
L T (W )* (W 0T )' ? L corresponding to a motif of order i + j. This motif corresponds to two 
chains of length % and j, respectively, emanating from the same root cell. To understand 
the impact of second order motifs we need to decompose this motif as illustrated in Figs. [7] 
and [8j While this is a challenging combinatorial problem, we show that the answer can be 
obtained by rearranging the terms in Eq. ( |23| . 

Each factor of the form L T (W°) ! (W 0T ) J L in Eq. d23b can be split by inserting I = H+@ 
between each occurrence of W° or W 0T , as in Eq. (24). The resulting expression can be 



used to identify the impact of motifs of order k on terms in the expansion of order i + j > k. 
The following Proposition, proved in Appendix [Ej formalizes these ideas. 

Proposition 4.1. Let H be the rank-1 orthogonal projection matrix generated by the unit 
N '-vector u ; H = uu T , © = I — H. For any N x N matrix K ; let 



Kr 



(Ker'K 



K0K 0K. 

V v ' 

n factors of K 



// the spectral radii ^(K) < 1 and $(K0) < I, then 

u T d-Kr 1 fi-K r r 1 u 



-i 



-i 



i - J2 uTk « u 1 + 2 uTK n© K - u 1 - uTk ™ u 



n=l 



n,m=l 



m=l 



(2f 



We will use Prop. 4.1 to derive a relation between second order motif strengths and 
mean covariances. Assuming that ^(Au'W ), $(IwW°0) < 1, and setting u = y/Nh and 
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K = Au/W°, applying Prop. 4.1 to Eq. ([7]) gives 

' oo > 

l-^A^l/W^L 



n,m=l 



(29) 



where 



W° 

n 



N 



m=l 



w°ew° • • • ew°, 

n-l 



n factors of W° 



w 



™> m fJn+m— 1 



w°ew° • • • ew° e w 0T ew 0T • • • ew 0T . 



n factors of W° 



factors of W 0T 



Keeping only terms in Eq. (23) which can be expressed as polynomials of second order 



in motif frequency and connection probability is equivalent to keeping only terms involving 



Wj (connection probability), W" (chain motifs) and Wjj (diverging motifs) in Eq. (29). 
This yields an expression which involves only first and second order motif frequencies: 



(C°°) 
C° 



1+ (nAiu) q dh 



N 



I - ( NAw) p - (NAw) 2 q ch 



(30) 



Fig. M shows that the approximation to the covariance given by Eq. (30) provides a 



significant improvement over the second order truncation approximation in Eq. (22) (See 



Fig. [6]). We emphasize that it requires only three scalars that summarize the statistics of 
the entire connection graph: the overall connection probability and the propensity of two 
second-order motifs. We offer a heuristic explanation for the effectiveness of the resumming 
theory based on spectral analysis of W in Appendix [HJ 
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Figure 9: A comparison of the mean correlation obtained using Eq. ^ (horizontal axes) 
to the resumming approximation in Eq. (30) (vertical axes). The diagonal line y — x is 
plotted for reference. Each panel corresponds to a different scaling of coupling strength for 
the same set of 512 adjacency matrices same as Fig. [6| which are sampled from SONET (see 



Sec. 3.2). The effective coupling strength is characterized by the Erdos-Renyi spectral radius, 
and is recorded at the top of each panel (see Sec. 2.4). Here, the open circle indicates the 
level of correlation expected from an Erdos-Renyi network with the same overall connection 
probability and strength. 



If we expand the denominator in Eq. (30 ) as a power series in 
and keep only terms which are linear in g<j iv , q^, we obtain 



; p + ^NA W y q ch 



(nAw) 



+ 



N{Awf 



(C°°) _ 1 
C° ' N(l - NAwp) 2 ' (1 - NAwp) 2 



:<?div + 



2N(Aw) 



1 - NAwp) 3 



q ch + h.o.t. 



(31) 



This linear relation was used to estimate the regression coefficients in Fig. [5] (b ar plot). 
Finally, we note that similar means may be used to derive versions of Eq. (|30[) that (unlike 



Eq. (31 )) retain a nonlinear dependence on motif frequencies, but keep motifs of either higher 
or lower order. To approximate the impact of motifs up to order r, we would keep terms 
with factors WjJ,Wj m where n, n + m < r in Eq. (29). For example, if we take r = 1, we 



estimate the mean covariance based only on the probability of occurrence of first order motifs 
- that is, connection probability. This is equivalent to estimating the covariance in idealized 
Erdos-Renyi networks where gdiv and g c h, and their higher order analogs are precisely zero. 



In Eq. (29), we set all terms involving W to zero save Wj, giving 



(C°°) 
C° 



1 



A ( 1 - NAwpY 



(32) 



This predicted mean correlation is indicated by the open dots in Fig. [9] (which accurately 
describes Erdos-Renyi networks as shown in Fig. [2]). This again demonstrates how motif 
structures exert a large influence over averaged spike correlations. 
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4.4 Correlations in external input 



As is natural for communication among different brain areas or layers Shadlen and Newsome 



f 998a , or as arises for certain sensory inputs Doiron et al. , 2004 , we next consider the case 



in which the network under study receives an additional noisy input that is correlated from 



cell-to-cell; in other words, the cells receive common input from upstream sources de la 



Rocha et al., 2007, Josic et al., 2009, Trong and Rieke, 2008 



total covariance structure 



'-NN 



We take this input to have 
•I), so that the variance of 



0) C = o\\ + a 2 x p input 
such external input (not to be confused with the implicitly modeled external noise £ (t) term 
in Eq. ([3])) to each neuron is fixed (ul] independent of p m P ut and the correlation coefficient 
of the inputs to all cell pairs is p m P ut Lindner et al. 2005, Marinazzo et al. , 2007| . Eq. ([7]) 



then has the form (see Trousdale et al. , 2012 for details, and an improvement applicable 



when the extra noise source is white): 

c°°=(i- Aw) 1 (c°i + i 2 C") (i - AW T 



= (i - aw) 1 [(c° + i 2 4) i + (i 2 4 P in p ut ) (i NN - 1)] (i - lw T ) 1 . 

In this case, the output variance in the absence of coupling (W = 0) predicted by the linear 
response theory is (C° + A 2 o~ 2 x ). Normalizing and multiplying by L T , L to arrive at an 



approximation of average correlation, and applying the ideas of Sec. 4^3, we find that 

(C°°) B 1 [l + (NAw) 2 q div ]+NB 2 



C° + A 2 a 2 



x 



N 



where 



By 



C° + A 2 a\{\ 



P 



i 2 ' 

(NAw)p - (NAw) 2 q ch 

input \ /)2^ r 2 ^input 

) n 71 a xP 



Bo 



C° + A 2 a 2 x C° + A 2 a 2 x 

(compare Eq. p0|)). We might ask how the presence of input correlations affects output 
correlations by calculating the change in output correlation resulting from input correlations 
of size p input : 



output 



(C°°) 



(C°°) 



C° + A 2 a 2 



x 

B 2 [N-1 



C° + A 2 a 2 x 
- (NAw) 2 q d 



pinput-o 



(33) 



N 



(NAw)p - (NAw) 2 q ch 



2 ■ 



For simplicity, we can linearize this expression in g^iv; Qch to examine the interaction between 
second order motifs and input correlations, giving 



Ap 



output 



?1 

N 



N — 1 



(NAw) 



[1 - {NAw)p} 2 [1 - {NAw)p} 



2(N-l)(NAw) 2 

> gdiv+ [\-{naw)pY qch 



(34) 
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If we assume that w ~ O (^) (so that the Erdos-Renyi spectral radius NAw ~ 0(1), see 
Sec. 2.4), then, asymptotically, gdiv only has an order O(j^) (negative) contribution, while 
g c h has a order 0(1) contribution to Ap output (note that the first term in (34) is also 0(1), 
which represents the "base" response to correlated input in a Erdos-Renyi network). This 
implies that in large networks, the chain motif is the most important motif in determining 
how input correlations will be transferred into network correlations - which will be "output" 
to the next area downstream. 



5 Impact of second order motifs in networks of excita- 
tory and inhibitory neurons 

Biological neuronal networks are composed of excitatory and inhibitory neurons (EI net- 
works). To treat this case, we next show how our theory extends to the case of networks 
composed of two interacting subpopulations. 



5.1 Results: linear dependence between mean correlation and mo- 
tif frequencies 

As in the previous section, we start by a numerical exploration of the contribution of motifs to 
network-averaged correlation in excitatory-inhibitory networks. We generated 512 networks 



using the degree distribution method, as described in Sec. |3.2| We evaluated network- 
averaged correlations using the full linear response expression given in Eq. ([7]). We then 
performed a linear regression analysis of the dependence of network correlations on motif 
frequency in networks of excitatory and inhibitory neurons, for the 20 second order motifs 
shown in Fig. |4| 

The linear regression gives a reasonable fit (R 2 > 0.69, see caption). Chain and diverging 
motifs contribute significantly to the network-averaged correlation, with some motifs having 
a higher effect than others. Moreover, mean correlations depend only weakly on converging 
motifs. Note also that these regression coefficients correspond to how we devise the axes in 
Fig. [3^3. There, we use these regression coefficients (or actually the theory predication for 



these values based on Eq. (44)) as weights to linearly combine multiple motif frequencies. 
Specifically, we take linear combinations of the 6 diverging, 6 converging and 8 chain motif 
frequencies (see Fig. |4j) respectively, giving the 3 axes in Fig. [3t 
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Figure 10: Linear regression coefficients between network-averaged correlations within differ- 
ent subpopulations (vertical axes) and motif frequencies (horizontal axes). The three panels 
correspond to averages across EE, EI and II cell pairs, from top to bottom. The linear re- 
gression coefficient between average correlation and 20 motif frequencies are represented by 
a dot (See Fig. [4] for the enumerated list of second order motifs in an EI network). Error bars 
represent a 95% confidence intervals. The vertical bars were obtained from the resumming 
theory described in Sec. 5^3 Here R 2 = 0.91, 0.86, 0.69 respectively. The spectral radius of 
AW under the Erdos-Renyi assumption is 0.25, and Ne = 80, Ni = 20, = 3.707 

so that pNeWe + pNjWi ~ giving approximate balance between average excitatory and 
inhibitory inputs. Parameters are given in Fig. |3j 
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5.2 First theory: second order truncation for network correlation 

To extend the theory developed in Sec. [4] to EI networks, we need to take into account 
distinct second order motifs. For instance, there are eight different types of three-cell chains 
(See Fig. Although this makes the notation more burdensome, the main ideas are the 
same. Indeed, following a similar approach, the theory can be extended to an arbitrary 
number of subpopulations. 

Consider a network of size N = Ne + Nj, consisting of Ne excitatory and iVj inhibitory 
neurons. Excitatory (resp. inhibitory) connections have weight We (resp. Wj), so that 
we > and Wi < 0. The connection probability from class X to class Y is pyx- Note that 
the connectivity matrix can now be written as 



W 



W EE W E i\ _ fw E W° EE WjW° EI 



W 7E Wjj) \w E W° IE wjW n 



where, Wjy and are respectively, the weighted and unweighted connection matrices 

between cells of class Y to cells of class X. The population sizes, weights, connection prob- 
abilities, and firing rates together determine the balance between excitatory and inhibitory 
inputs that cells receive (see Discussion). Throughout this section we use networks with 



same parameters and architecture as in Fig. 10 



First, define the N x 2 block-averaging matrix L by 

L E \ = [1 NB)1 /N E 

where Ijvm is the N x M matrix of ones. We define the two-population analogs of the 
orthogonal projection matrices H and © by 

H = LD 2 L T , and © = I — H, (35) 

where D 2 is the 2x2 matrix 

D - ( Ne 
2 "V0 N I/ 

Direct matrix multiplication shows that if X is a matrix with block form 



X 



X^£ X.EI 

X/_e X77 



where Xy^ is an Ny x Nz matrix, then L T XL is a 2 x 2 matrix of block-wise averages of 
X, that is, 



<X) fl — L XL — l t^ ieLe W l - I {Xje) (X//) 



We will make use of the empirical average connection strength matrix M given by 

M = L T WL = ( WePee WiPei ) 
\wePie wipn J 
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To examine the dependence of mean correlations within a block on the frequency of second 
order motifs, we consider the block-wise average of the covariance matrix C°°, 



(C°°)b = L T C°°L 



'P$b) (Cf/> 
(CTe) (Cfj) 



Excitatory and inhibitory connection weights need not be equal. It is therefore necessary 
to consider motif frequencies simultaneously with connection weights. For instance, the 
contributions of a length two chain passing through an excitatory or inhibitory intermediary 
cell, such as motifs 13 and 14 in Fig. |4| are not necessarily equal and opposite in sign. Their 
contributions are dependent on the ratio we/wj. To account for this, we define motif strength 
matrices Qdiv, Qcon, and Q c h as follows. 

The strength of diverging motifs expected in an Erdos-Renyi network is given by 

Qd£ = MD 2 M T = (L T WL)D 2 (L r W T L) 
= L T WHW T L 

_ / N e (wePee) 2 + N^wtPei) 2 N e w e PeePie + NiWjPeiPiA 
\N e w 2 e peePie + NjwjpEipji N e (wePie) 2 + N 1 {w I p II ) 2 ) ' 

Here multiplication by D 2 converts average individual motif strength (e.g. w e PeePie) to 
average total motif strength. The empirical average of the strength of diverging motifs is 
given by 

Qtotal = L T WW T L 

(WeeWIe) + (Wb/WIj) (WeeWJe) + (Ws/Wj,) 
(W IE Wl E ) + (Wj/WL) (W ie WJe) + (Wj/WJj) 

Hence, we can write the expected total strength of diverging motifs in excess of that expected 
in an Erdos-Renyi network as 

n — I ^div Qdil 1 def n total n ER 
^div — I qIE Qll J — ^div ^div 

KNEwlq^ + Njwjq^ 1 N E w%q^f + Njwjq 1 ^ J (36) 
= L' WW' L - L T WHW T L 
= L T W(H + 0)W T L - L T WHW T L 
= L T W0W T L. 

XY,Z 1 /-tuT-O ihtOT \ 
Qdiv = Jj- ( W XZ^Yz) ~ PXZPYZ 

represents the probability of observing a diverging motif to cells of classes X, Y from a cell 
of class Z in excess of that expected in an Erdos-Renyi network. 



where 
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It is important to note that the matrix Qdiv contains motif strengths in excess of what 
would be expected in an Erdos-Renyi network (i.e., number of occurrences, scaled by con- 
nection weights and probability of occurrence), while the scalars gaiv still correspond to 
probabilities. 

The strengths and frequencies of converging motifs can be expressed similarly, giving 



Q c 



f'W'WL - M T D 2 M = l/W T 0WL 



qEE qEI\ 
Vcon Vcon I 

IE O u 1 

Vcon Vcon/ 

/ N E w 2 E q E f E + Njwlqgf 1 N E w E w iq E ^ E + iWurf/ 
\N E w EWl q^ E + NjWEWjq^ 1 N E wjq^ E + Njwjq 1 ^ 



(37) 



where 



XY,Z 
^con 



1 

N~ z 



wf z w° YZ ) 



PZXPZY, 



represents the probability of observing a converging motif to cells of class Z from cells of 
classes X, Y in excess of that expected in an Erdos-Renyi network. 
Finally, for chain motifs, 



Q 



ch 



Qch Qch 



Q IE 



ch 



Q" 



ch 



def 



L 1 W Z L - MD 2 M = L T W0WL 



( N E w 2 E q EEE + NjWEWjqg 11 ^ N E w E Wjq^ EE + Njwjq 1 ^ 
V N E w 2 E q EEI + NjWEWjqg 11 NeWeW^ 1 + Njwjqg 1 



(38) 



where 



ZYX 
Qch 



iVv 



(W ZY W YX ) 



PzyPyx 



represents the probability of observing a length two chain motif beginning at a cell of type 
X and terminating at a cell of type Z, passing through a cell of type Y, in excess of what 
would be expected in an Erdos-Renyi network. 

A truncation of Eq. (16) at second order gives an initial approximation of the block-wise 
average (C(0))b- 



(C°°) B /C { 



L + h.o.t. 



= L T 1 + A (W + W T ) + A 2 (W 2 + W 2T + WW T ) 
w [L T L + A (M + M T ) + A 2 (MD 2 M + M T D 2 M T + MD 2 M T ) 
+ A 2 (Qch + Qch + Qdiv)] • 



(39) 



Terms not involving a second order motif matrix correspond to the contributions of motifs up 
to second order in a purely Erdos-Renyi network — for example, MD 2 M gives the expected 
contribution of length two chains in an Erdos-Renyi network. 

Fig. 11 compares the second order truncation given by Eq. (39 ) with the mean correlations 
obtained from the entire series in Eq. (16). The correlations in the network of inhibitory 
and excitatory cells can be appreciable, even when the spectral radius of matrix AW is 
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much smaller than one. However, the second order truncation of Eq. (39) gives a poor 
approximation of network-averaged correlations. 



Excitatory-Excitatory Average Excitatory-Inhibitory Average Inhibitory-Inhibitory Average 

Approx. Approx. Approx. 




0.05 



0.2 All motifs 




-0.04 



0.05 All motifs 




-0.04 -0.02 All motifs 



Figure 11: Block- wise average correlations obtained from the second order truncation of 



Eq. (]39|). On the horizontal axis are the average correlation between cells from the given 

and the approximate value obtained from the truncation is 
x corresponds to perfect agreement between 



classes (calculated from Eq. (|7J) 
given on the vertical axis. The diagonal line y 
the true value and the approximation. The spectral radius of AW under the Erdos-Renyi 



assumption is 0.33 (see Sec. 2.4). Other network parameters are the same as in Fig. 10 12 
and 13 Here R 2 = 0.33, 0.02, 0.0005 respectively for the three panels, confirming that the 
truncation approach gives a poor prediction of network correlation. Here, the open circle 
indicates the level of correlation expected from an Erdos-Renyi network with the same 
overall connection probabilities and strengths. 



5.3 Improved theory: resumming to approximate higher-order 
contributions to network correlation 

As in the case of a single population, we can improve our prediction of mean correlations by 
accounting for the contributions of second order motifs to all orders in connection strength. 
The equivalent of Eq. ( 23 ) has the form 



i9^h = i-i/C^L = i i+i L T W i (W T )^L, (40) 

C ^ i,j=0 

where W and L are as defined in the previous section. 

First, we generalize Prop. 4.1 to the case of two populations (for a full proof, see Ap- 
pendix [F]) 

Proposition 5.1. Let H = UU T be a orthogonal projection matrix generated by a N x M 
matrix U, whose columns are orthonormal vectors. Define = I — H. For any N x N 
matrix K ; define K ra as 

K„ = K0K 0K. 

V 

n factors of K 
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If spectral radius ^(K), $(K0) < 1, we have 



U J 'fI - KrVi - K 



u 



i - J2 uTK n u 1 + E U T K n 0K^U I - £ U T K^U 



n=l 



n,m=l 



m=l 



We now apply this Proposition to the expression in Eq. (40). Let U 



LD 



1/2 



(41) 



and 



U T XU = D* /2 L T XLD* /Z for any matrix X, so that H has the form given in Eg. (J35J). In 
addition, let K = AW, and assume that $(iW), *(i6W) < 1. Then Prop. 

(C 00 )^/^ = L T (I - iW)- x (I - iW T )- ! L 



lV2 



5.1 



gives 



-i 



n,m=l 



I - ^i n L T W n LD 2 ( D 2 X + ^ i" +m L T W n 0W^L 

n=l / 
' oo > 

I - £ v4 m D 2 L T W^L 



(42) 



m=l 



where 



w n = wew ew . 

n s ^ , 

n factors of W 



As in the single population case, we can discard all terms in Eq. (42) which do not correspond 
to second order motif frequencies. This means that in the first and third set of brackets in 
Eq. (42 ), we discard any terms containing W n with n > 3, and for the middle set of brackets, 
we discard terms containing W n @W^ for n + m > 3. This gives 



(C^b/C « (i - iMD 2 - i 2 Q ch D 2 ) 1 (d 2 1 + i 2 Q div ) (i - AD 2 M T - A 2 B 2 Q, 

' (43) 

Figure 12 illustrates that this approximation is a great improvement over that given by 
truncating the expansion at second order (compare with Fig. 11). Again, we note that 
our approximation requires only knowledge of the overall connection probabilities among 
excitatory and inhibitory cells, and the frequency of second-order motifs. 
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Figure 12: Predicting block- wise average correlations from resumming theory. The 
horizontal-axis is the average correlation (in a certain block) from original linear response 
expression of covariance matrix (Eq. [7]); the vertical axis the quantity from resumming the- 
ory (Eq. (43)). The diagonal line y = x is plotted for r eference. The spectral radius of 



AW under Erdos-Renyi assumption is 0.33 
the same as in Fig. 10 Here R 2 



(see Sec. 



2.4). Other network parameters are 



0.93, 0.88, 0.87 respectively for the three panels. The open 
circle indicates the level of correlation expected from an Erdos-Renyi network with the same 
overall connection probability and strength. 



Expanding the inverses in Eq. (43) in a power series, we can again obtain an approxima- 



tion to block average correlations to linear order in Q c h, and Q 



div; 



(CHb/C « ( I — v4MD 2 



+ A 2 Q 



D, 1 



A z Q ch I - AB 2 M 



I - AM T D 



-l _ 



A 2 Ql 



div 



I - AD2M 1 



(44) 



As in the single population case, each entry of the 2x2 matrix on the right hand side of 



Eq. (44) gives an approximation to the block-averaged correlations expressed in terms of the 
scalars q^ x 9di v ' > P rov iding an analytical estimate to the regression coefficients plotted in 



Fig. 10 An example with uniform connection probability is given in Appendix G We note 



that in general all sub-types of diverging and converging motifs affect each block-wise average 



correlation. This is not predicted by the second order truncation in Eq. (39). Somewhat 
counterintuitively, q^ E and q^ 1 (index 3 and 4 in Fig. [4J can contribute negatively to 
(Cee) as shown in Fig. 
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As in the single population case, we can also retain contributions to mean correlation 
which can be expressed as (nonlinear) functions of first order motifs (connection probabilities) 
only. This follows from setting all Q terms in Eq. (43) to zero, yielding the following 



approximation for mean correlation in the two-population analog of Erdos-Renyi networks: 



(C°°) B /(7 = (I - AMD 2 



-1 



D, 1 



I - iD 9 M T 



-1 



(45) 
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This predicted mean correlation for an Erdos-Renyi network is shown by an open dot in 



Fig. [12] the deviations from this value illustrate that motif structures can both significantly 
increase and decarese average network-correlations (see also Fig. [2] and [9]). 



6 Heterogeneous networks 

For simplicity in the previous sections we assumed a homogeneous network of neurons com- 
posed of cells with identical firing rates, power spectra and response properties. As a result, 
the diagonal matrices C° and A in our key expression for correlations, Eq. ([6]), were scalar 
matrices, C°I and Al, and could be factored - leaving the connectivity structure of the net- 
work to determine the correlation value. In particular, this structure impacts correlations 
via the matrix products of adjacency matrices in Eq. (16). 

Biological neural networks are heterogeneous. Even if the neurons are identical before 
they are coupled, any heterogeneities in the coupling structure will lead to different neurons 
firing with different rates and hence power spectra (cf. Zhao et al. 2011 , which we return 



to in the Discussion). Moreover, as a consequence, they will also have different levels of 
responsivity. In sum, the matrices C° and A will not have identical entries on the diagonal. 



Consequently, the equivalent of the expansion Eq. (16) takes the form 



(46) 



i,j=0 



As a result, a purely graph theoretic interpretation of the terms in the expansion - that is, 
one based on the connectivity alone - is no longer possible. Recall that as in Eqns. (|9]|TIJ, 
motif frequencies are associated with terms like W t (W T ) J . Here the powers (AW)* are of 
a "weighted" connection matrix (with weights corresponding to the responsivity of different 
cells). Moreover, C° is no longer a scalar matrix and in general does not commute with W, 
introducing additional complications that we address below. 

In this section we discuss how to extend our results to the case of heterogeneous neural 
populations. 



6.1 Performance of the homogeneous approximation 

A first attempt at coping with heterogeneity is to hope that it is unimportant, and to apply 
the homogeneous theory of Eq. (16) naively. To do this, we need to choose approximate 



(scalar) values for the power spectrum C° and responsivity A that we will apply to all of 
the cells, by plugging into the homogeneous network formula given by Eq. (16). We choose 



the C° as the unadjusted power spectrum (due to the normalization ([2]), the actual value 
of C° is not important), rather than using one homogeneous value for all cells. For A, we 
use the geometric mean of the A^ which are found from self consistent equations similar to 
Eq. ([5]) (see Trousdale et al. 2012 for details). The results of this naive application of the 



homogenous theory are shown in panel A of Fig. 13 Although general trends are captured, 
this approach does not give an accurate approximation. 
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We note that here coupling is relatively strong (Erdos-Renyi formula for spectral radius 



giving 0.4, see Sec. 2.4). With weaker coupling (spectral radius=0.3), and hence less network- 
driven heterogeneity, the homogeneity assumption provides improved approximations in the 
heterogeneous case (R 2 measure 0.78, 0.67, 0.67 for EE, EI and II average correlation re- 
spectively, data not shown). To quantify the heterogeneity of cellular dynamics across the 
networks, we compute the coefficient of variation CV of A\ and C® averaged over network 
samples (CV = (standard dev./mean)). For spectral radius=0.4, the CV is 0.35 and 0.38 
for Ai and Cf, respectively, while at spectral radius=0.3, the CV is 0.23 and 0.26 for Ai and 
Cf. At such levels of heterogeneity, we clearly need a more systematic approach, and we 
develop this next. 

6.2 Heterogeneous theory 



We are prevented from applying Proposition 4.1 to Eq. (46) in the heterogeneous case, due 
to the presence of the factor C° in the middle of the terms on the right hand side. To deal 
with this difficulty we use the substitution C° = (C°) 1 / 2 (C ) 1 / 2 , which is possible because 
the power spectrum is non-nonnegative. We can then rewrite Eq. ^ as 

c(o) « (c ) 1/2 (c°r 1/2 (i - Awr^cy/^c ) 172 ^ - w^Xr^c )- 1 / 2 ^ ) 1 / 2 

= (C ) 1 / 2 (i - (C°)- 1 / 2 AW(C ) 1 / 2 )" 1 (i - (Cy^W^C )- 1 / 2 ) -1 (C ) 1 / 2 , 

(47) 

where we are again evaluating all quantities at u = 0, so that F = I. 

Let Ao be the geometric mean of Ai, which we choose to normalize responsivity so that 
weighted quantities will have the same units (see below). We can then define an effective or 
"functional" connection matrix, which has the same units as W : 

W = (C )- 1/2 AW(C°) 1/2 /i , (48) 



so that Eq. (47) becomes 

(C°) 1/2 (I - ioW)- 1 ^ - I W T )- 1 (C°) 1/2 . (49) 

This expression will be much easier to study, as the diagonal matrix C° no longer appears in 
the middle. We can thus expand the two terms, (I — v4 W) _1 and (I — A W T )~ 1 to obtain 



an expression that has a form analogous to that of Eq. (16). 

The only difference in the present case of a heterogeneous network is that the definition 
of motif frequency and connection probability will involve weighted averages. For example, 
the entries of W are scaled version of entries of the connection matrix W: 



A n X /W; 
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As an example, a diverging motif i ^ — A; — >• j has a weighted contribution of the form 



AiAj 



A 2 



=±=W ik W jk . 



(50) 



The ratio (AiAj) / Aq in Eq. (50) quantifies the relative responsivity of the recipient cells. 
Hence, a particular diverging motif will be weighted more strongly (and provide a greater 
contribution to average correlation) if the recipient cells are more responsive to inputs. 
Similarly, C°, C®, C° corresponds to the variance of spike counts in long time windows for 
the uncoupled cells: the weight is determined by the variance of the projecting cell divided 
by the geometric mean of the variance of the recipient cells. These observations agree well 
with intuition: more responsive cells will be more strongly correlated by a common input, 
and "source cells" with larger total variance (C®) will lead to a diverging motif with larger 
impact. 



The motif frequencies gdiv and g c h, are defined by Eq. (18) and (20) upon substituting W° 
with W° = W/w in the single population case. In the case of two populations, the matrices 
Qdiv and Q c h are defined by Eq. (36) and (38), with W replaced by W. These weighted 



motif frequencies could be estimated experimentally in an active network, via recordings 
from neurons known to participate in three cell motifs. The advantage of our resumming 
theory remains that everything needed is based only on the statistics of a relatively small 
number of cell motifs, rather than higher order information about the connectivity graph. 
Due to the normalization in Eq. (|2]), the two outside diagonal matrix factors in Eq. 



(49) cancel in the definition of p, and we can use the following matrix to calculate corre- 



lation coefficients (noting that we continue to approximate Cu by Cf when performing the 
normalization in the denominator of Eq. (|2l) 



(I - A W) _1 (I - A W 



T\-l 



(51) 



Eq. (51) is the heterogeneous analog of the expression for the correlation C^/C stud- 
ied in Sections lil and [B] (see Eq. rt7|), including diagonal contributions. Applying the motif 



resumming theory exactly as in these sections yields corresponding approximations of mean 
correlation for these systems. In particular, for the two population case, the mean correlation 



Pee-i Pei i Pn w * n ^ e estimated based on Eq. (43) with the substitution of weighted motif 



counts Qdiv and Q c h (note that one must still adjust the approximation to account for diago- 
nal terms — see Appendix [A}. In panel B of Fig. 13 , we plot this prediction of p#J, p^f, pff 
compared with the full expression for correlation in the heterogeneous case (via Eqns. Q 
and (J6|). Accounting for dynamical heterogeneity across the network by defining weighted 
second order motifs provides a reasonably accurate prediction of average correlation. 
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Figure 13: A Mean correlation from homogeneous resumming theory (horizontal axis) 
comparing with that from full linear response theory (Eq. (J7J)) (vertical axis). B Mean 
correlation from heterogeneous resumming theory (horizontal-axis) comparing with that from 
full linear response theory (vertical-axis). The diagonal line y = x is plotted for r eference. 
The spectral radius of AW under Erdos-Renyi assumption is 0.40 (see Sec. 2.4). Other 



network parameters are given in the caption of Fig. 10 The coefficients of determination, 
R 2 , are 0.56, 0.44, 0.35 in panel A and 0.88, 0.81, 0.73 for panel B. 



7 Comparisons with IF simulations 



In the preceding sections, we have shown how network-wide correlation can be approximated 
based on coupling probability together with the frequency of second-order motifs. In assess- 
ing the accuracy of this motif-based theory, we have compared the predictions of the theory 
with the value of correlation given Eq. (|6]), which is exact in the sense that it precisely 
includes the contribution of motifs at all orders and therefore gives an exact description of 



linearly interacting point process models Pernice et al. , 2011, Hawkes, 1971a 



When using our theory to predict the impact of motifs on mean correlation for networks 
of IF neurons, we are making an additional approximation in describing integrate-and-fire 



dynamics with linear response theory Trousdale et al. , 2012, Ostojic et al. 2009, Lindner 



et al. , 2005 . We now directly test the performance of our motif-based predictions in IF 



networks - thus probing how errors at each of the two layers of approximation combine. 



Specifically, in Fig. 14 we compare the block-wise mean correlations from IF simulation and 
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the predications obtained using Eq. (43). These simulations are for the same networks used 



m Fig. [3j We find that our theory gives a good quantitative match with LIF simulations, 
despite the multiple approximations involved. Thus, our theory predicts trends in the impact 
of motif frequencies on network-wide correlation. 



Excitatory-Excitatory Average 



Excitatory-Inhibitory Average 



Inhibitory-Inhibitory Average 



Simulations 



Simulations 



Simulations 



0.1 



0.05 



0.02 



' ./-V 



0.1 0.2 Nonlinear 

Resumming Theory 



0.05 



0.1 Nonlinear 
Resumming Theory 




0.02 Nonlinear 

Resumming Theory 



Figure 14: Average correlation from integrate and fire neuron simulations compared with 
the predictions of Eq. (43). The vertical axis is the mean correlation coefficient calculated 



from simulations of the same 265 excitatory-inhibitory networks studied in Fig. [3j The 
horizontal axis is the prediction for mean correlation from the resumming theory based on 
empirical connection probability and second order motif frequencies. The diagonal line y = x 
is plotted for reference. The spectral radius of AW under the Erdos-Renyi assumption is 



0.33 (see Sec. 2.4). Coefficients of determination R 2 are 0.91, 0.87, 0.79 respectively for the 
three panels. 



8 Discussion 



Summary: Predicting network-wide correlation from three cell mo- 
tifs 

We studied the impact of the graphical structure of neural networks - characterized by con- 
nectivity motifs - on their dynamical structure, characterized by correlated activity between 
cells. As shown in Fig. |2j varying the frequency of such motifs can strongly impact corre- 



lation, over and above the overall level of connectivity in a network. Following Zhao et al. 



2011 , we focus on the three types of motifs that involve two connections each: the diverging, 



converging and chain motifs. 

We chose a standard spiking neuron model, the integrate-and-fire (IF) neuron, in con- 
structing our recurrent networks (Sec. [2|. For IF neurons Lindner et al. 2005, Trousdale 



et al., 2012 , LNP neurons Beck et al. , 2011 



and other neuron models such as linearly inter- 
acting point process model Hawkes, 1971a|b , one can apply linear response approximations 
or, in the case of Hawkes models, solve exact equations - see "Neuron models" below) to get 
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an explicit expression for the pairwise correlations depending explicitly on the connectivity 
matrix (see Eq. Q). 

We expand this expression in a series, where each term has clear correspondence to certain 



graphical structures Pernice et al. 2011 , Rangan 2009a . In particular, second order terms 



correspond to second order motifs. Importantly, we show that contributions of higher order 
terms can be estimated using the frequencies of second order motifs along with the connection 
probability. 

For systems with correlations well-approximated by Eq. ^ and assuming homogeneous 
cellular dynamics, we find that the frequency of converging motifs, over and above that 
expected in Erdos-Renyi networks, will have no effect on mean correlation in systems (if the 
connection probability and other motif frequencies are fixed) . Meanwhile diverging and chain 
motifs will contribute positively. In networks of excitatory and inhibitory neurons, the three 
types of motifs are subdivided according to the type of the constituent neurons. However, 
average correlations between cells of a given type are still given in terms of diverging and 
chain motifs (see Fig. [3]). 

We first analyzed networks in which we made a strong homogeneity assumption on the 
dynamical properties of the uncoupled neurons. The resumming theory we develop approx- 
imates the contributions of higher order motifs in terms of the frequency of second order 
motifs. In Sec. |6j we extended our theory to heterogeneous networks. In such cases, the con- 
tribution of one instance of a certain motif to the total motif frequency will be additionally 
weighted by the relative responsiveness of the neurons composing the motif which receive 
input, as well as the baseline variance of the cells. Overall, our results can be regarded as a 
general estimator of mean correlation given motif statistics up to second order. 

To test our theory numerically, we generated random networks of IF neurons with fixed 
statistical (i.e., expected) connection probability, but different second order connectivity 
statistics (see Sec. [3]). Simulations show that the theory provides an accurate description 
of network correlations (see Fig. [I4| ), despite the additional error introduced by the linear 
response approximation of activity (Eq. Q). We also compared the resumming theory to 
the direct evaluation of the linear response theory (Eq. ([6])) which takes account of the full 



graphical structure. The close match between the two (Fig. |9j Fig. 12) shows that second 
order motifs capture much of the dependence of mean correlation on network connectivity. 
Moreover, Figs. [2j |9j 12 demonstrate that variations in the frequency of second order mo- 
tifs produce changes in network correlation that dominate those expected from the small 
variation in empirical (i.e., realized) connection probability from one network to the next. 

Beyond the resumming theory, we also considered two other ways of simplifying the ex- 
pansion of Eq. Q. The first was a truncation in connection strength at second-order (in 



powers of A and W; Sees. 4.2, 5.2). This eliminates all contributions due to motifs of length 
two or less. Except in very weakly connected networks, this is a poor approximation: Al- 
though the contributions of higher order motifs decay exponentially in interaction strength 
Aw, their number also grows exponentially with Np. Thus, motifs of all orders could, in 
principle, contribute substantially to average correlations. Our second truncation was the 
Erdos-Renyi approximation of mean correlation. This yielded predictions for mean correla- 
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tion that included contributions from motifs of all orders. These predictions therefore depend 
only on connection probability (See Eqns. (32, 45)), and should be valid when there is very 
little structure in the connectivity graph compared to an idealized Erdos-Renyi network. 



Neuron models 



Our motif-based theory can be applied to a variety of different neuron models. The starting 
point of our theory, the expression for pairwise cross-spectra given by Eq. ([6]), arises in a 
number of settings. 

The main tool in our approach is linear response theory. This connects our methods to 
IF models (see Sec. 2.3). Importantly, Eq. ( |6j) also ar ises as an exact expression for linearly 



interacting point process, or Hawkes models Hawkes, 1971b, Pernice et al.[ |2011| [2012 
this case, 

C(u) = (I - i(w)W)~ 1 Y(I - i(u;)W T )~ 1 



In 



where A{uS) is the Fourier transform of the interaction filters and Y is a diagonal matrix of 
firing rates. Therefore our analysis can be directly and exactly applied to this setting. We 
note that the Hawkes process is a linear-Poisson (LP) model. Such models are commonly 



applied in theoretical neuroscience Dayan and Abbot, 2001 . 



Relationship to other studies on motifs and networks correlation 



Our methodology is very similar to the previous study of Zhao et al. 2011 . They considered 



the impact of diverging, converging and chain motifs (along with the reciprocal connection 
motif) on the ability of an excitatory recurrent network to stay synchronized. They found 
that prevalence of the converging motif decreases synchrony, prevalence of the chain motif 
increases synchrony, and that the diverging motif has no significant effect on synchrony. The 
difference between their results and ours can be understand from the different dynamical 
regimes considered. Zhao et al. 2011 considered perturbations from two extreme cases: 
perfect synchrony, and evenly distributed asynchronous oscillators. In contrast, we studied 
the asynchronous regime but allowed the activity of the cells to be correlated. Hence, our 
methods also differ: we use a linear response approach valid for strong internal noise, and 



2011 performed a linear stability analysis 



weak to intermediate coupling, while Zhao et al. 
of coupled oscillator equations. 

Many other studies have also examined the relationship between graphical features and 
spike correlations in networks. Roxin 2011 studied how spike-time correlation changes 



when one interpolates the degree distribution of a network from the binomial distribution 
(corresponding to Erdos-Renyi networks) to a truncated power law distribution. He found 
that increasing the variance of the out-degree distribution will increase the cross-correlation, 



which is consistent with our results for the diverging motif (see Eq. (12)). 



Pernice et al. 2011 and Trousdale et al. 2012 



tures on pairwise correlation for a number of network classes. Pernice et al. 2011 



studied the influence of connection struc- 

used the 



Hawkes neuron model, which as noted above, leads to very similar expressions for correlation 
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as those studied here and by Trousdale et al. 2012 . Both approaches relate pairwise cor- 
relations to certain graphical structures of increasing order (i.e., motif size). In particular, 



Pernice et al. 2011 obtained an expression for average network correlation in terms of the 



mean input and common input (Eq. (22) in Pernice et al. , 2011] ) for regular networks with 



uniform connection probability, while Trousdale et al. 2012] (Eq. (25) in Trousdale et al. 



2012| ) considered correlations in networks where only in-degree was fixed. Both are special 
cases of our resumming theory (Note that according to Eq. (12) and (13), a fixed in (resp. 
out) degree is equivalent to q COQ = (resp. gdiv = 0) and g c h = 0). 

A major contribution of the present study is to show that effects of "higher order" 
graphical interactions (i.e., motifs including more than three cells) can be approximated in 
terms of the frequency of second order motifs and the overall connection probability. This 
allows a systematic treatment of network-averaged correlation for a broader range of network 
connectivities. 



Limitations 



When applied to integrate and fire neuron networks, our our analysis relies on the validity 
of the linear response approximation. In Section [7j we demonstrated this validity for a 
particular firing regime, for the class of random networks studied here. For more on this 
Trousdale et al. 2012 and Ostojic and Brunei, 120111. We note that one avoids 



issue, see 



this issue entirely when considering Hawkes processes (see above). 

We also assumed that the spectral radius of the total coupling matrix is less than one, 
in order to expand Eq. ^ into a series, for which each term can be attributed to a different 
motif. From this point, our methods rely on our ability to predict the impact of higher 
order motifs on network correlation based on the frequency of second order motifs. We 
demonstrated that our resumming theory can successfully make this prediction for classes of 
networks generated in two ways: via two-sided power law distributions, and via the SONET 
method (see Sec. 3.2). However, for certain connectivity matrices, our resumming theory can 



produce large errors. An example pointed out to us by Chris Hoffman [personal communica- 
tion] is W° containing a ^JpN x y/pN block of 1 entries and but taking value everywhere 
else; here, p is the overall connection probability. Note that this matrix corresponds to a 
graph with one fully connected group and another fully isolated group. Such disconnected 
architecture and strong inhomogeneity may be features that produce large spectral radii 
$(W8), and hence large errors when applying our theory (se e also Appendix [5] a nd for a 
study on the dynamics of inhomogeneous clustered networks, |Watts and Strogatz , 1998] ) . 

Three other factors limit the generality of our results. First, in order to separate the 
effect of network structure from that of cellular dynamics we initially assumed homogeneous 
firing rates and identical neurons. For many real neural networks, this may be a poor 
approximation. Thus, in our analysis of heterogeneous networks, we calculate weighted 
motif strengths, assuming "a priori" knowledge of the (heterogeneous) cellular firing rates 
and responsivities. A full theory would rather begin by predicting this heterogeneity based 
on network properties. Intriguingly, Zhao et al. 2011 show that certain motif frequencies 



can be a powerful source of such heterogeneity, a potentially important fact that we neglect. 
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Second, for simplicity we only considered long-time window spike count correlation. How- 
ever, our analysis can be easily generalized to study correlation at any time scale, and for 
any u in Eq. ([6]). A third, and final limitation of our analysis is that we only studied pairwise 
correlation coefficients. This may not adequately describe the network dynamics or reveal 
the full importance of certain motifs. 



Extensions and connections with neurobiology and computation 



Intriguingly, experiments have shown that motif statistics in actual neural networks deviate 
from the Erdos-Renyi network model, opening a possible role of multicellular motifs. For 
example, Song et al. 2005 found such deviations in connectivity between excitatory cells 



in visual cortex. Sporns and Kotter 2004 also suggest that increased frequencies of certain 
motifs may be a common feature of neural circuits. Moreover, Haeusler et al. , 2009 studied 
the different motif statistics from two experimental data of connectivity structures of laminae 
in cortex, showing deviations from expectations under the Erdos-Renyi model without the 
laminar structure (implying as a possible origin of non-trivial motif frequencies). Our study 
can be applied to suggest a possible consequence of these experimental findings for levels 
of spike-time correlation (in the latter case generalizing it to apply to multi-group networks 
corresponding to the different laminae). 



Renart et al. , 2010 studied pairwise correlations in balanced excitatory- inhibitory Erdos- 
Renyi networks, and found that balanced neuronal networks will have average correlations 
that decrease rapidly to zero with network size . Our results suggest that, with increased 
propensity of certain second order motifs, EI networks can potentially produce significantly 
larger correlations value (see Fig. [2]). However, to apply our findings to the case considered 
by Renart et al. 2010 1 , we would need to make several extensions, including allowing full 



heterogeneity arising from recurrent connections and checking that linear response theory 
works in their strongly coupled, balanced dynamical states. In our setting, "balanced" can 
be defined as that the total mean input received for each cell is 0, exactly or on average 

This needs to be compared 



Raj an and Abbott 




2006 


in Renart et al. 




2010 . 



For future studies, we also hope to develop a theory that predicts not just the mean 
correlation strength across a network, but its variance across cell pairs. This variance has 



been studie using theoretical and experimental approaches Renart et al. , 2010, Ecker et al. 



2011 , and it would be interesting to describe how it depends on connection motifs. We will 
also try to predict the heterogeneity in cellular dynamics caused by motif frequencies, as 



referred to above - incorporating, for example, the result emphasized in |Zhao et al. , 2011 



that variability in input to different neurons depends on converging motifs (Eq. ( 13 )). Finally. 



we note the important connections between correlations and coding efficie ncy jZohary et al. 



1994 Abbott and Dayan] 1999 Salinas and Sejnowski 2000 Josic et al. , 2009 Schneidman 



et al. 2006 (see also Introduction). The recent work of Beck et al. 2011 set up a direct 
connection between the graphical structure of networks and their coding efficiency, using a 
similar linear response calculation of covariance matrix for linear nonlinear Poisson (LNP) 



neurons. The present results cold be used with the approach of Beck et al. 2011 to further 
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link the statistics of motifs to properties of signal transmission in neural networks. 



A Approximating p avg from (C°°) 

The average covariance across the network, (C°°) can be used to approximate p avg . Here 
we describe two such approximations. First, when the uncoupled neurons have equal firing 
rates, and the perturbation from recurrent coupling is weak, the diagonal terms C°° will 
be close to the unperturbed values, C°. In this case, subtracting diagonal terms from the 
average, we have that 

^-f^--V— (52) 
P \&> Nj N-l 

In a second, more accurate approximation, we assume permutation symmetry between 
neurons within one population. Also, in our networks, self-connections are allowed and occur 
with the same probability as other connections. These will lead to identical (marginal) 
distributions for each entry in the covariance matrix C°°, excepting the diagonal entries 
which are shifted by a constant of C° due to each neuron's own unperturbed variance (this 
corresponds to the term proportional to I if one were to expand Eq. ([7j as a series — see 



Eq. (16)). This suggests that C°° has the form 

C°°nC l + cl NN , 

where Inn is the N x N matrix of all ones, and c is a constant. If this holds, then, using 
diagonal entries of this matrix as normalization, we obtain 

C„ + (C~> - C„IN 



The two approximations given in Eqns. ( |52|[53 ) are approximately equal for small cor- 



relations. We will use Eq. (52) to exhibit the linear dependence between mean correlation 



coefficient and motif frequency (such as linear weights in Fig. |3] Part B) and Eq. (53) for 
quantitative predications (all numerical plots). 

For networks consisting of an excitatory and inhibitory population, we have analogs of 
Eqns. ( 52p3 ) for each of the 4 population blocks of the covariance matrix, 



and 



Pxx 



avg „ (CJ? X ) 1 \ N X 

Pxx \ C° N X J Nx-V 

p™y « <C*y(0)>/C(0), (54) 

(Cxx) - C /N x 
C + {C XX )-C Q /Nx 

(' 

(Co + (Cxx) - C /N X )(C + (C YY ) - Co/Ny) 



PXY ~ / ( 55 ) 
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where X e {E,I}. 



B Proof of bound on gdiv, Qch for one population 



We here prove the inequalities in Eq. (15). Noting that we may write 

gdiv = E e [w; fc wy - E e [wy e £ [w° j , 

(with similar expressions holding for q con , q ch ) it is sufficient to show E e [W° fc W° fc ] < p. Note 

HW°(W°) T H = (iV 2 E e [W° fc W° fc ] ) H, (56) 



where H is defined in Eq. (17). Let ||-|| F be the Frobenius norm, which is sub- multiplicative. 
We then have 

||HW°(W°) T H|| F < ||HW || F ||(W°) T H|| F < ||H|| F ||W || F ||(W°) T || F ||H|| F . (57) 



Since ||H|| F = 1, ||W°|| F = y Ylii j^Wi,j) 2 = Ny/p, the above inequality, together with (56), 

gives E e [Wg fe Wj fc ] < p. With the fact that E e [W°J E e [W°J = p 2 , we have the bound in 
Eq. (pi). 



In the second inequality of (57), we have used ||HW || F < ||H|| F ||W || F . The neces- 
sary and sufficient condition for achieving equality is the equality condition in the following 
Cauchy inequalities for all i, j 



E^wj,) s(e*&) (e«/) 

< k J \ k / \ k J 



Equality holds exactly when W®j = for all k, I - that is, each column of W° has 
same values (either all ones or all zeroes). The equality condition for the first inequality in 
Eq. (57) is identical. To generate an example graph that achieves equality for a specified 
overall connection p, we simply set a fraction p of the columns to be all ones, and set the 
remaining columns to zero. 

Similarly, for converging motifs we can show E^W^W^-] < p with the same equality 
condition as above, and for chain motifs, E e [W° fc W° J < p. However, the equalities for 
<?div, Qcon cannot hold simultaneously, and the equality for q ca cannot be achieved. 



C Graph generation methods 

Here we present more details on how we generated network samples with fixed connection 
probability, but different frequencies of second order motifs. We used two methods. 
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First, the degree distribution method consists of initially generating a sample of in and 



out degrees from a truncated power law distribution with density, following Zhao et al. 



2011 




< d < L x 
Li<d<L 2 
otherwise 



(58) 



where d is the in or out degree (see also the configuration model Roxin, 2011, Newman 

2001] ). The two marginal distributions of in and out degree 



2003, Newman and Watts 



are then coupled using a Gaussian copula with correlation coefficient p to generate the in 
and out degree lists. The parameters p, Li/L 2 ,L 2 , 71 > 0,72 < are independently and 
uniformly sampled for each network, separately for in and out degrees (their ranges are listed 
111 Table Q. d and C 2 are chosen so that f(d) is continuous at L\ and the mean of the 
degree distribution is normalized to Np (the same value for both in and out degrees), where 
p the fixed connection probability across network samples and N is the network size. 

We then use the degree lists to calculate a probability for each possible connection from 
cell j to cell i, which is proportional to G?f<i° ut . All these iV 2 probabilities are scaled so that 
the resulting average for the total number of connections is the same as the quantity iV 2 p stat . 
Here p stat is the target connection probability that we aim to achieve in these samples; we 
recall that the "empirical" connection probability achieved in a given graph is notated as p. 





P 


7i 


72 




L 2 /N 


min 


-1 


0.25 


-2.25 





0.7 


max 


1 


2.25 


-0.25 


1 


1 



Table 1: Range of parameters of the truncated power law degree distribution. 



For the excitatory-inhibitory case, we generate four degree lists d^, rf^ ut , df, dj ut , again 
according to the marginal distributions (58). We therefore need a four dimensional Gaussian 
copula. Again, parameters for the power law distribution and the correlation coefficient 
matrix of the copula are randomly chosen. Using these degree list, we can then generate 
each of the four blocks (defined by cell type) of adjacency matrix in the same way as in the 
single population case. This method allows us to sample from the whole extent of motif 
parameters (3 in single population case, 20 in two population case). 

For the single population networks studied in Sec. [4j we generated additional network 
samples via the SONET method Zhao et al. 2011 . The idea of the algorithm is similar to 



maximum entropy models. Given only the connection probability and second order motif fre- 
quency, we try to generate the most "random" network satisfying these moment constraints. 
However, instead of using the Gibbs distribution for the connections as in maximum entropy 
method, we use a dichotomized (N 2 dimensional) Gaussian distribution. We then generate 
samples of excitatory only networks from the complete range of possible motif frequencies, 
as described in Eqns. (14)-(15), using the SONET algorithm. 
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Network samples generated using both methods cover the range of motif frequency ob- 
served experimentally in cortical circuits Song et al. 2005, Zhao et al. , 2011 as shown by 
Table [2] Here, we list motif frequency values as (gdiv> <Zcon> <Zch)/(j>(l ~ P)) f° r excitatory only 
networks, and {q EE,E , qc Q E,E , q E ^ E ) I {Pee{^- — Pee)) in excitatory-inhibitory networks (since 
2005| only recorded from excitatory neurons, see definition at Eq. (36)) . 



Song et al. 







SONET: E only 


degree: E only 


degree: EI 




experiment 


min 


max 


min 


max 


min 


max 


diverging 


0.033 


0.001 


0.913 


0.015 


0.181 


0.018 


0.295 


converging 


0.044 


0.001 


0.838 


0.016 


0.189 


0.022 


0.279 


chain 


0.022 


-0.192 


0.248 


-0.068 


0.095 


-0.082 


0.110 



Table 2: Range of motif frequencies in network samples. "E only" ( or "EI") means excitatory 
only networks (or excitatory- inhibitory networks). "SONET" and "degree" refer to the two 
methods of generating networks. 



D Compensating for fluctuations in empirical connec- 
tion probability 

In order to isolate the impact of higher order motifs, all network samples should have the 
same empirical connection probability p, i.e. the same number of connections, as we now 
explain. The algorithms used to generate networks produced samples with slightly different 
empirical connection probabilities. These fluctuations impact the linear relationship between 
motif frequencies and average correlation, as p factors into the regression coefficients (See 



Eq. (31)). When attempting to determine the regression coefficients from data, fluctuations 
in p affect the linear trend. 

To address this issue, we can scale motif frequencies by weighting them so that the 
contribution of a motif to mean correlation does not have a regression coefficient depending 
on p. For example, suppose that the theory predicts a linear relation of the form 

(C°°) 

- = /o(p) + fch{p)q c h + fdiv(p)qdiv 
Then we may define auxiliary motif frequencies 

q , = hip) 

fx{Pstat) 

and also replace fo{p) with fo(p s t a t) so that the theoretically-predicted regression relationship 
becomes 

£, Q = fo(Pstat) + fch(Pstat)q c h + fdiv {Pstat) ?div • 
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Here, we have adjusted the definition of a motif frequency in order to account for finite- 
size fluctuations in the empirical connection probability p. The linear regression fit to the 
quantities q' cii) q'^ v is much improved, achieving an R 2 measure of 0.99, up from 0.8 in the 



case where we did not account for such fluctuations. In Fig. [15j we show the scatter plots 
exploring the relationship between motifs and mean correlation after performing this scaling, 
and observe the same key trends of mean correlation as in the unadjusted Fig [5j strong and 
positive dependence on chain and diverging motifs, and minimal dependence on converging 
motif. These trends are now presented even more clearly. 
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Figure 15: Reproduction of Fig. [5] where we have scaled motifs to account for fluctuations 
in p, as described in the text. 



E Proof of Proposition 4.1 



We will make use of the following lemma, which may be verified by direct computation. 

Lemma E.l. Let {x„}„>i, {y m } m >i, {z>nm\n,m>\ be sequences which converge absolutely 
when summed, and also satisfy 



Then, 



E 

71=1 



<1, 



oo 
m=l 



< 1. 



OO / OO 



E E IK)=E(E 

i=l (ni,...,n fe )e{j} V=l 



i=l \ra=l 



E E 

M = 1 (n 1 ,...,n k+ i)e{i] 
{m 1 ,...,m l+1 )&{j} 



n 

vs=l 



X n s I z n k+1 m l+ 



i Y[Vmt 



oo / oo 

EE 

i=0 \n=l 



E 

vn,m=l 



oo / oo 
j=0 \m=l 
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where the sum over (n 1; . . . , n^) G {i} denotes a sum over all ordered partitions of i, 
(ni, . . . , nk), of length 1 < k < i, with each n\ > 1. 

A general result is that, for any matrix A, there exists a sub-multiplicative matrix norm 
(||XY|| < ||X||||Y||) arbitrarily approaching the spectral radius A), that is ||A|| < \I/(A)-|- 



Horn and Johnson, 1990 . Since $(K0) < 1, we can choose a sub-multiplicative matrix 



norm ||-||_x satisfying ||K0||^ < 1. As all matrix norms are equivalent, there also exists a 
constant c satisfying ||-|| 2 < c||-||a Horn and Johnson, 1990 . For n > 1, we have 



|u T K„u| = |u T (K0 



,n-l 



Ku 



< llul 



2 • 



(K0 



,n-li 



■ IlKlUllul 



< clKKer-^UIlKlla 



< cIlKllallK© 



m-i 



Thus Y^=i u K n u converges absolutely, as it is bounded above in absolute value by a 
convergent geometric series. Since u T K„u = u T K^u, this also implies the convergence of 
Ylm=i uT Km u - Absolute convergence of Y^ m =i uT K n 0K^u is proved similarly by noting 
that for n,m > 1, 



u 1 K n 0K> 



= u T (K0) n K T (0K T ) u 

< ||u|| 2 ||(K0r|| 2 ||K T || 2 ||(K0; 
<c 2 ||K|| 2 ||fK0) n llAll (K0) m - ] 



2 U 2 



< c 2 ||K|U||K0 



|A| 
\n+m—l 



Now, suppose that |$^°=i u T K n u| < 1. We will remove this assumption shortly. We can 
now expand the right hand side of Eq. (28) as 



oo / oo 



i=0 \n=l 



1 + u T K,0Kju 

i,j>l 



oo / oo 



j=0 \m=l 



(59) 



Using Lemma E.l, we have that 



(59) 



oo / oo 



+ 



£ £u T K„u 

\n=l , 

oo / oo 

£ £" T K„u 



OO / oo 



i=0 \n=l 



j=0 \n=l 

oo 

u T Ki©Kju 



\i,3>l 



oo / oo 



E Ea:» 

3=0 \n=l 



!+e e ri uTK »." 

i=l (ni,...,n fc )e{i} \s=l 



! + E E IK K ». U 

3=1 (n u ...,m)e{j} \t=i 



(60) 



M = l (ni,...,n fe+ i)e{i} 
(mi,...,m ;+ i)e{i} 



vt=l 



47 



Next, note that the product term []j =1 u T K„ s u, where (n 1; ...,n fc ) G {i}, is acquired 
by distributing across sums H + in u [K (H + 0)f Ku and taking the unique term in 
this where factors of H = uu T divide the i factors of K in to k blocks K ns of size n s joined 
by factors of 0. Therefore, summing over all possible ordered partitions, we have 



u T K*u = u T (K(H + 9)) w Ku = II uTR ^ u ) 

(ni,...,n fc )e{i} \a=l 

Similarly, 

u T K*0(K T ) J 'u = u T (K(H + 0)) < - 1 K0K T ((H + &)K T y- l u 



e n uTr - u u T K nfc+1 0K^ i+i u n uTk ^ u 



(ni,...,n fc+1 )e{«} \a=l 
(mi,...,m l+1 )e{j} 



vt=l 



Using Eqns. (60 62), we have that 



(61) 



(62) 



© = E uTlCu E uT ( KT ) Ju + J2 u T K*0(K T ) J 'u 



,i=0 
oo 



0=0 



i,j>l 



E u T K J H(K T )^u+ E u T K i 0(K T )^u 



»J>0 



»J>0 



= E uTKi ( kT ) J u usin § I = H + 

i,i>0 

= u T (I-K)- 1 (I-K T )- 1 u, using tf(K) < 1, 

where we have used u T 0(K T ) J u = u T K'0u = on the second line of the equation. 

Lastly, we will eliminate the assumption X^Li uT K n u| < 1, which was used to establish 
). To see that this can be done, let z be a complex number, and replace K by zK. in 
), giving 



28 



Eq. 
Eq. (128 



u T a-zK)- i a-zK T )- 1 u 



-i 



-i 



1 - £ z n u T K n u 1 + X) ^ +m u T K n 0K^u 1 - z m u T Klu 



n=l 



n,m=l 



m=l 



(63) 



For sufficiently small < S < 1, \z\ < 8 we have that 



X z n u T K n u 



< 5^2\u T K n u\ < 1. 



n=l 



n=l 
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The absolute convergence of the series J2^=i uT K„u, Y^=i uT K n u and Y^ m =i uT K„0K^u 
then guarantees that Eq. (63) holds on \z\ < 5, and furthermore, that the right hand side of 
(63) is an analytic function of z on \z\ < 5. 



Finally, note that the left hand side of (63) is an analytic function of z on \z\ < 1/\I/(K), 
since the matrix inverses may be expanded as power series on this range. Since the two sides 
are equal and analytic on \z\ < 5, the right hand side must also be defined and analytic on 
\z\ < 1/\&(K), and equal the left hand side on this range. Since 1/\&(K) > 1, we may take 
z = 1 in Eq. (63), recovering Eq. (28). 



F Proof of Proposition 5.1 



The proof of Proposition |5. 1| is identical to the one dimensional case Prop. 4.1 except that 
we need an entrywise argument to eliminate the assumption ||$^Li U T K n U||2 < 1 that 
enables the expansion of the right hand side of Eq. (41). Let z be a complex number, and 



replace K by zK in Eq. (41), giving 



U T (I-*K)- 1 (I-zK r )- 1 U 

f oo \ - 

I - ^"U T K„U ( I + ^ +m U T K„9KlU ) ( I - ^ m U T K^U 



-i 



n=l 



n,m=l 



m=l 



-1 



(64) 



Note that the series Y^=i U T K n U absolutely converges in 2-norm. This follows from writing 



|U T K n U|| 2 < ||U T || 2 ||(K0) n_1 || 2 ||K|| 2 ||U|| 2 < cllU' 



, 2 ,|K0||riK|| 2 ||U|| 2 , 



and use the assumption that \& (K0) < 1 and one associated sub-multiplicative norm 
|| KG) || a < 1. c is a constant such that ||-|| 2 < c||-||a- For sufficiently small < 5 < 1, 
\/\z\ < 5, we have that 



|^ n U T K n U|| 2 < <5^||U T K n U|| 2 < 1. 



n=l 



n=l 



With this condition, we can establish identity (64) in the same way as in the one dimensional 



case Prop. 4.1 for |z| < S. Furthermore, it is similar as above to show all of the matrix series 



on the right hand side of Eq. ( 64 ) are absolute converge in 2-norm and therefore entry-wise 



absolute converge. Hence, we have that each entry of the right hand side of Eq. (64) is an 
analytic function of z on |z| < 5. 



On the other hand, the left hand side of Eq. (64) is entry- wise analytic in z (each entry is 



indeed a rational function of z) on \z\ < 1/^(K). Since the two sides are equal and analytic 
on |z| < 5, the right hand side must also be defined and analytic on \z\ < 1/\I/(K). Note 
that 1/\&(K) > 1, so we may set z — 1, yielding Eq. (41). 
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G Expression of the linear dependence between (Cee) 
and motifs 

For excitatory-inhibitory networks, the (linearized) relationship between motif frequencies 
and averaged correlations is given by Eq. (|44j). Here, we evaluate the terms in this expression 
in a special case, to yield an explicit expression for the linear dependence of correlation on 
motif frequencies. We will express the block-wise average covariances (Cee) m terms of the 
20 individual second order motifs. Along with the (approximate) definition of correlation in 



Eq. (52), this gives the linear weights. We note that this is how we obtain the specific linear 
weights used in Fig. [3]B. 

For simplicity, assume the special case that all connection probabilities are identical 
(Pxy = P for all X, Y e {E,I}). Define the weight of all excitatory (resp. inhibitory) 
connections into a cell as 

\i E = pN E w E (resp. Hi = pN^i) 
and the new weight of all connections into a cell as 

A* = Ve + Vi- 
lli addition, define rj to be the strength of total common input to a cell pair in an Erdos-Renyi 
network: 

r] = N E w 2 E p 2 + N lW jp 2 . 



Noting that we have 
it is simple to show that 



MD 2 ) 2 = (MD 2 ) 



Similarly 



-i A 

I - AMB 2 ) = I + ^ (MD 2 ) • 

1 — Afi 



-i A 
I-AB 2 M) =1 + — (D 2 M) . 

i — An 
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Then, from Eq. (44) we get the linear dependence on motifs (treat p as fixed) as 



(C EE )/C° ~ i 2 



1 - Ap T 

l-Ap 



n EE i ~ A 2 Api(l — A/ij) £J ~ 2 



~ 2 l-i M/ (l-w+^f 



+ 2A 



1 - Ap 



Q 



EE 
ch 



+ 2A 2 
+ 2i 
+ 2i 2 



~ 2 1 - A/i/ (1 - Ap E )Api + (1 - Apj)Ap E 



El 



I- An 

2 Ani ( 



i-W + iV 2 f 

Vch 



<3 



Si 
ch 



;i - Apy 



Apj (1 - Ap E )Api + (1 - Api)Ap E 



Mr 



l-Ap (1 - A//) 2 

where Q^,Q^ were defined in Eqns. ( 36p8 ) as 



Vch> 



Vdiv 



N EW E q div + N lWl q div , 



Q 



XY 
ch 



N E w x w E q ch 



XEY 



Api 



1 - Ap 



Q 



11 

div 



+ N!Wxw iq ^ Y . 



(65) 



Since each of the quantities , are clearly linear in second order motif frequencies, 
Eq. (65) gives a linear relation between second order motif frequencies and block- wise av- 
eraged correlation in the two population network. Eq. (65) is the two population analog of 
Eq. (31) in the single population case. 



H Intuition for why the resumming approach can pro- 
duce accurate results 



From Eq. (29) and Eq. (42) , we see that the error of resumming theory is determined by the 
tail that we dropped in two series. With respect to the coupling strength order of magnitude 
(w or w E ,wj), these are geometric series. Therefore the sum of the tail series is controlled 
by the leading term, that is 

i 3 L T wewewL, i 3 L T wew T ewL, i 3 L T wewew T L, 



This shows that the resumming theory has third order accuracy in the effective interaction 
strength Aw. 

Another important factor affecting the accuracy of the resumming theory is the spectral 
radius \I/(@W@) = \&(W@) = ^(©W) (equalities follow from writing W under basis of 
projections and H), which is related to how fast the terms in the series converge to 0. 
There is a simple intuition of \I/(W@) for Erdos-Renyi networks. For single population 
networks, note that (asymptotically, for large N) W°@ = W° — W°H is a matrix with 
i.i.d. entries of mean 0. According to the Circular Law, all eigenvalues will asymptotically 
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uniformly distributed within a circle about the origin. Comparing with the spectra of W° 
in Sec. 



2.4 



Abbott, 2006 



multiplying by effectively removes the single dominant eigenvalue Raj an and 



W°e as 



2.4). Such 



The removal of this dominant eigenvalue will reduce the spectral radius o 
compared to W° by a factor of y/N (w\Jp{l — p)N compared to pNw, see Sec. 
reductions of $(W0) approximately occur in single population networks and at the blocks 
of W© in two population networks, even though those networks are non-Erdos-Renyi. This 
intuition may help understand why resumming theory works much better than truncation 
theory: the tails of series in W@ may be much lighter than those in W. 



Imag 4 



-e- 

20 
Real 



Figure 16: The spectra of single population Erdos-Renyi network, the larger circle and small 
circle on the right are expected locations of the bulk spectra and single eigenvalue calculated 
from the asymptotic formula in text above and Sec. 2.4 The excitatory only network has 
100 neurons, and p = 0.2, w = 1. 
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Symbol Description 



Vi,Ti, Elj, Oi Membrane potential, membrane time constant, leak reversal potential, 

and noise intensity of cell %. 

Ei,<Ji Mean and standard deviation of the background noise for cell i. 

Vth, v r , r re f Membrane potential threshold, reset, and absolute refractory period for 

cells. 

i/)(v),Vt, A t Spike generating current, soft threshold and spike shape parameters for 



the IF model Fourcaud-Trocme et al. 2003 



fi(t),r]i(t) Synaptic input from other cells in the network, and external input to 

cell i. 

T s,i, T D,i Synaptic time constant and delay for outputs of cell i. 

yi(t) Spike train of cell i. 

Wy The j — y i synaptic weight, proportional to the area under a single 

post-synaptic current for current-based synapses. 

3ij(t) The j — > i synaptic kernel - equals the product of the synaptic weight 

and the synaptic filter for outputs of cell j. 

Cjj(r) The cross-correlation function between cells i,j defined by CV,-(r) = 

cov(yi(t + r), Uj(t)). 

C = C(0) The cross-spectrum matrix evaluated at frequency. Unless noted 

otherwise, all spectral quantities are evaluated at frequency. 

N y .(t,t + r), Pij(r) Spike count for cell i, and spike count correlation coefficient for cells 
i, j over windows of length r. 

Ti, Aj(t),C°j Stationary rate, linear response kernel and uncoupled auto-correlation 

function for cell i]. 

Kjj(t) The j — > i interaction kernel - describes how the firing activity of cell 

i is perturbed by an input spike from cell j. It is defined by K^-(t) = 

(A*Jii)(*)- 

y"(t), Cfj(t) The n th order approximation of the activity of cell i in a network which 

accounts for directed paths through the network graph up to length n 
ending at cell i, and the cross-correlation between the n th order approx- 
imations of the activity of cells 

g(t),g(u) 9^) is the Fourier transform of g(t) with the convention 



/oo 
e-^g{t)dt 
-oo 



E e [-] Empirical average, ^ ' or w J2ij,k ' depending on context 

R 2 coefficient of determination, i.e. the square of correlation coefficient p 2 

W° Adjacency matrix 

H H = jjInn 

e © = i — h 53 



Table 3: Notation used in the text. 
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